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