DECLARE FUNCTION chi2# (X AS DOUBLE, df AS INTEGER)
DECLARE FUNCTION poz# (z AS DOUBLE)
DECLARE FUNCTION xinvchi2# (P AS DOUBLE, df AS INTEGER)

'**************************************************************************
'*                  CHI2 AND INVERSE CHI2 FUNCTIONS                       *
'* ---------------------------------------------------------------------- *
'*  Reference:                                                            *
'*                                                                        *
'*  http://www.fourmilab.ch/rpkp/experiments/analysis/chiCalc.html        *
'*                                                                        *
'* The following JavaScript functions for calculating chi-square probabi- *
'* lities and critical values were adapted by John Walker from C imple-   *
'* mentations written by Gary Perlman of Wang Institute, Tyngsboro.       *
'* Both the original C code and the JavaScript edition are in the public  *
'* domain.                                                                *
'* ---------------------------------------------------------------------- *
'* SAMPLE RUN:                                                            *
'*                                                                        *
'*    CHI2 Law:                                                           *
'*    ========                                                            *
'*  Value of variable: 100                                                *
'*  Degree of freedom: 100                                                *
'*                                                                        *
'*  Probability =  .4811916845279567                                      *
'*                                                                        *
'*  Verify Inverse Law:                                                   *
'*    X =  100.000000148592                                               *
'*                                                                        *
'*                                  Basic Version By J-P Moreau, Paris.   *
'*                                           (www.jpmoreau.fr)            *
'**************************************************************************
DEFDBL A-H, O-Z
DEFINT I-N

CLS
PRINT
INPUT " Value of variable: ", X
INPUT " Degree of freedom: ", idf
PRINT
P = chi2(X, idf)
PRINT " Probability = "; P
PRINT
PRINT " Verify:"
Y = xinvchi2(P, idf)
PRINT "   X = "; Y

END 'of main program


FUNCTION chi2 (X AS DOUBLE, df AS INTEGER)
' Adapted From:
'Hill, I. D. and Pike, M. C. Algorithm 299
'Collected Algorithms for the CACM 1967 p. 243
'Updated for rounding errors based on remark in
'ACM TOMS June 1985, page 185

DIM even AS INTEGER  'parity of df
DIM True AS INTEGER
DIM False AS INTEGER
DIM lnpi AS DOUBLE
DIM ipi AS DOUBLE
    PI = 3.14159265359#
    bigx = 200#
    True = 1: False = 0
    lnpi = LOG(PI ^ .5)
    ipi = 1# / LOG(PI)
    IF X <= 0# OR df < 1 THEN
      chi2 = 1#
      EXIT FUNCTION
    END IF
    a = .5 * X
    IF INT(df / 2) = df / 2 THEN
      even = True
    ELSE
      even = False
    END IF
    IF df > 1 THEN
        Y = EXP(-a)
    END IF
    IF even = True THEN
      s = Y
    ELSE
      s = 2 * poz(-(X ^ .5))
    END IF
    IF df > 2 THEN
        X = .5 * (df - 1)
        IF even = True THEN
            z = 1#
        ELSE
            z = .5#
        END IF
        IF (a > bigx) THEN
            IF even = True THEN
                e = 0#
            ELSE
                e = lnpi
            END IF
            c = LOG(a)
            DO WHILE (z <= X)
                e = LOG(z) + e
                s = s + EXP(c * z - a - e)
                z = z + 1
            LOOP
            chi2 = s
        ELSE
            IF even = True THEN
                e = 1#
            ELSE
                e = ipi / (a ^ .5)
            END IF
            c = 0#
            DO WHILE (z <= X)
                e = e * (a / z)
                c = c + e
                z = z + 1
            LOOP
            chi2 = c * Y + s
        END IF
    ELSE
        chi2 = s
    END IF

END FUNCTION

FUNCTION poz (z AS DOUBLE)
'----------------------------------------------------------------------
'POZ  --  probability of normal z value

'Adapted from a polynomial approximation in:
'Ibbetson D, Algorithm 209
'Collected Algorithms of the CACM 1963 p. 616

'Note:
'This routine has six digit accuracy, so it is only useful for absolute
'z values < 6.  For z values >=  6.0, poz() returns 0.0.
'----------------------------------------------------------------------
DIM Y, X, w AS DOUBLE
     zmax = 6#
    IF z = 0# THEN
      X = 0#
    ELSE
      Y = .5 * ABS(z)
      IF Y >= zmax * .5 THEN
        X = 1#
      ELSEIF Y < 1# THEN
        w = Y * Y
        X = ((((((((.000124818987# * w - .001075204047#) * w + .005198775019#) * w - .019198292004#) * w + .059054035642#) * w - .151968751364#) * w + .319152932694#) * w - .5319230073#) * w + .7978845605929999#) * Y * 2
      ELSE
        Y = Y - 2
        X = (((((((((((((-.000045255659# * Y + .00015252929#) * Y - .000019538132#) * Y - 6.769049860000001D-04) * Y + .001390604284#) * Y - .00079462082#) * Y - .002034254874#) * Y + .006549791214#) * Y - .010557625006#) * Y + .011630447319#) * Y - _
 9.279453341000001D-03) * Y + .005353579108#) * Y - .002141268741#) * Y + .000535310849#) * Y + .999936657524#
      END IF
    END IF
    IF z > 0# THEN
        poz = ((X + 1) * .5)
    ELSE
        poz = ((1 - X) * .5)
    END IF
    bigx = 20
    
END FUNCTION

FUNCTION xinvchi2 (P AS DOUBLE, idf AS INTEGER)

    epsilon = .000001
    chimax = 100000#
    xminchisq = 0#
    xmaxchisq = chimax
    IF P <= 0# THEN
      xinvchi2 = xmaxchisq
      EXIT FUNCTION
    ELSE
      IF P >= 1# THEN
        xinvchi2 = 0#
        EXIT FUNCTION
      END IF
    END IF
    chisqval = (1# * idf) / (P ^ .5)
    WHILE ABS(xmaxchisq - xminchisq) > epsilon
        chisave = chisqval 'chisqval is modified by call to chi()
        IF (chi2(chisqval, idf) < P) THEN
            xmaxchisq = chisave
        ELSE
            xminchisq = chisave
        END IF
        chisqval = (xmaxchisq + xminchisq) * .5
    WEND
    xinvchi2 = chisqval

END FUNCTION