'***************************************************
'*          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