'*************************************************** '* Statistical distributions * '* ----------------------------------------------- * '* This program allows computing several dis- * '* tributions: * '* 1. binomial distribution * '* 2. Poisson distribution * '* 3. normal distribution * '* 4. normal distribution (2 variables) * '* 5. chi-square distribution * '* 6. Student T distribution * '* ----------------------------------------------- * '* REFERENCE: "Mathematiques et statistiques By H. * '* Haut, PSI Editions, France, 1981" * '* [BIBLI 13]. * '* ----------------------------------------------- * '* SAMPLE RUN: * '* * '* Statistical distributions * '* * '* Tutorial * '* * '* 1. Input distribution number:" * '* * '* 1: binomial distribution" * '* 2: Poisson distribution" * '* 3: normal distribution" * '* 4: normal distribution (2 variables)" * '* 5: chi-square distribution * '* 6: Student T distribution" * '* * '* 2. Define the parameters of chosen distribution * '* * '* 3. Input value of random variable" * '* * '* * '* Input distribution number (1 to 6): 3 * '* * '* Normal distribution * '* * '* MU=mean * '* S =standard deviation * '* X =value of random variable * '* * '* MU = 2 * '* S = 3 * '* X = 5 * '* * '* Probability of random variable = X: .0806569 * '* Probability of random variable <= X: .8413447 * '* Probability of random variable >= X: .1586552 * '* * '*************************************************** DEFINT I-N DEFDBL A-H, O-Z DIM B(4) 'used by normal distribution CLS PRINT " Statistical distributions" PRINT PRINT " Tutorial" PRINT PRINT " 1. Input distribution number:" PRINT PRINT " 1: binomial distribution" PRINT " 2: Poisson distribution" PRINT " 3: normal distribution" PRINT " 4: normal distribution (2 variables)" PRINT " 5: chi-square distribution" PRINT " 6: Student T distribution" PRINT PRINT " 2. Define the parameters of chosen distribution" PRINT PRINT " 3. Input value of random variable" PRINT PRINT INPUT " Input distribution number (1 to 6): ", ich CLS IF ich < 1 OR ich > 6 THEN PRINT " Error: Invalid choice!" END END IF 'call appropriate subroutine ON ich GOTO 100, 200, 300, 400, 500, 600 100 PRINT PRINT " Binomial distribution:" PRINT PRINT " P=probability of success" PRINT " N=number of trials" PRINT " X=value of random variable" PRINT INPUT " P = ", p INPUT " N = ", n INPUT " X = ", x GOTO 700 200 PRINT PRINT " Poisson distribution" PRINT PRINT " MU=mean" PRINT " X =value of random variable" PRINT INPUT " MU = ", xm INPUT " X = ", x GOTO 700 300 PRINT PRINT " Normal distribution" PRINT PRINT " MU=mean" PRINT " S =standard deviation" PRINT " X =value of random variable" PRINT INPUT " MU = ", xm INPUT " S = ", s INPUT " X = ", x GOTO 700 400 PRINT PRINT " Normal distribution (2 variables)" PRINT PRINT " MX=mean of X" PRINT " MY=mean of Y" PRINT " SX=standard deviation of X" PRINT " SY=standard deviation of Y" PRINT " RO=correlation coefficient" PRINT " X =value of 1st random variable" PRINT " Y =value of 2nd random variable" PRINT INPUT " MX = ", xm INPUT " SX = ", sx INPUT " MY = ", ym INPUT " SY = ", sy INPUT " RO = ", ro INPUT " X = ", x INPUT " Y = ", y GOTO 700 500 PRINT PRINT " chi-square distribution" PRINT PRINT " NU=number of degrees of freedom" PRINT " X =value of random variable" PRINT INPUT " NU = ", nu INPUT " X = ", x GOTO 700 600 PRINT PRINT " Student's T distribution" PRINT PRINT " NU=number of degrees of freedom" PRINT " X =value of random variable" PRINT INPUT " NU = ", nu INPUT " X = ", x 700 ON ich GOSUB 1000, 2000, 3000, 4000, 5000, 6000 'print results PRINT a\$ = " Probability of random variable = X: " B\$ = " Probability of random variable <= X: " c\$ = " Probability of random variable >= X: " d\$ = " Probability of random variables = X,Y: " IF ich = 4 THEN PRINT d\$; fxy: END IF ich <> 6 THEN GOTO 800 PRINT " Prob(-X<=random variable<=X) = "; bt PRINT B\$; px PRINT c\$; qx END 800 PRINT a\$; fx PRINT B\$; px PRINT c\$; qx END 1000 'Binomial distribution subroutine '**************************************************** '* INPUTS: * '* p: probability of success * '* n: number of trials * '* x: value of random variable * '* OUTPUTS: * '* fx: probability of random variable = X * '* px: probability of random variable <= X * '* qx: probability of random variable >= X * '**************************************************** q = 1# - p: t = p / q n1 = n + 1 fx = q ^ n'fx=prob(0) px = fx IF x = 0 THEN GOTO 1200 FOR i = 1 TO x fx = (n1 - i) * t * fx / i px = px + fx NEXT i 1200 qx = 1# + fx - px RETURN 2000 'Poisson distribution subroutine '**************************************************** '* INPUTS: * '* xm: mean * '* x: value of random variable * '* OUTPUTS: * '* fx: probability of random variable = X * '* px: probability of random variable <= X * '* qx: probability of random variable >= X * '**************************************************** fx = EXP(-xm) px = fx IF x = 0 THEN GOTO 2100 FOR i = 1 TO x fx = fx * xm / i px = px + fx NEXT i 2100 qx = 1# + fx - px RETURN 3000 'normal distribution '**************************************************** '* INPUTS: * '* xm: mean * '* s: standard deviation * '* x: value of random variable * '* OUTPUTS: * '* fx: probability of random variable = X * '* px: probability of random variable <= X * '* qx: probability of random variable >= X * '**************************************************** 'define coefficients B(i) B(0) = 1.330274429# B(1) = -1.821255978# B(2) = 1.781477937# B(3) = -.356563782# B(4) = .31938153# pp = .2316419# y = (x - xm) / s'reduced variable zy = EXP(-y * y / 2#) / 2.506628275# fx = zy / s 'calculate qx by Horner's method t = 1# / (1# + pp * y) po = 0 FOR i = 0 TO 4 po = po * t + B(i) NEXT i po = po * t qx = zy * po px = 1# - qx RETURN 4000 'normal distribution (2 variables) '**************************************************** '* INPUTS: * '* xm: mean of x * '* ym: mean of y * '* sx: standard deviation of x * '* sy: standard deviation of y * '* ro: correlation coefficient * '* x: value of 1st random variable * '* y: value of 2nd random variable * '* OUTPUTS: * '* fxy: probability of random variables * '* = (X,Y) * '**************************************************** s = (x - xm) / sx'1st reduced variable t = (y - ym) / sy'2nd reduced variable r = 2# * (1# - ro * ro) vt = s * s + t * t - 2# * ro * s * t vt = EXP(-vt / r) s = sx * sy * SQR(2# * r) * 3.1415926535# fxy = vt / s RETURN 5000 'Chisqr subroutine '**************************************************** '* INPUTS: * '* nu: number of degrees of freedom * '* x: value of random variable * '* OUTPUTS: * '* fx: probability of random variable = X * '* px: probability of random variable <= X * '* qx: probability of random variable >= X * '* gn: Gamma(nu/2) * '**************************************************** 'Calculate gn=Gamma(nu/2) xn2 = nu / 2 gn = 1# n = INT(xn2) IF nu = 2 * n THEN GOTO 5200 'Calculate gn for nu odd IF nu = 1 THEN GOTO 5100 FOR i = 1 TO 2 * n - 1 STEP 2 gn = gn * i / 2# NEXT i 5100 gn = gn * 1.7724528509# GOTO 5300 5200 'Calculate gn for nu even IF nu = 2 THEN GOTO 5300 FOR i = 1 TO n - 1 gn = gn * i NEXT i 5300 'Calculate fx fx = x ^ (xn2 - 1) aa = EXP(-x / 2#) bb = gn * (2# ^ xn2) fx = fx * aa / bb 'Calculate px fa = (x / 2#) ^ xn2 fa = fa * aa / (xn2 * gn) fm = 1#: px = 1#: i = 2 5400 fm = fm * x / (nu + i) px = px + fm IF fm * fa > 1E-08 THEN i = i + 2 GOTO 5400 END IF 'Accuracy is obtained px = px * fa qx = 1# - px RETURN 6000 'Student's T distribution '**************************************************** '* INPUTS: * '* nu: number of degrees of freedom * '* x: value of random variable * '* OUTPUTS: * '* bt: probability of random variable * '* between -X and +X * '* px: probability of random variable <= X * '* qx: probability of random variable >= X * '**************************************************** c = .6366197724# x = x / SQR(nu) t = ATN(x) 't=theta IF nu = 1 THEN bt = t * c GOTO 6500 END IF nt = 2# * INT(nu / 2) IF nt = nu THEN GOTO 6200 'calculate A(X/NU) for nu odd bt = COS(t) tt = bt * bt vt = bt IF nu = 3 THEN GOTO 6100 FOR i = 2 TO nu - 3 STEP 2 vt = vt * i * tt / (i + 1) bt = bt + vt NEXT i 6100 bt = bt * SIN(t) bt = (bt + t) * c GOTO 6500 'calculate A(X/NU) for nu even 6200 tt = COS(t) ^ 2 bt = 1#: vt = 1# IF nu = 2 THEN GOTO 6300 FOR i = 1 TO nu - 3 STEP 2 vt = vt * i * tt / (i + 1) bt = bt + vt NEXT i 6300 bt = bt * SIN(t) 6500 px = (1# + bt) / 2# qx = 1# - px RETURN 'End of file distri.bas