#! /bin/sh
echo ########
echo # editorial remark, 11 Nov 1995, ehg@research.bell-labs.com
echo # See also http://netlib.bell-labs.com/netlib/toms/322.Z.
echo # Matthew Belmonte <mkb4@Cornell.edu> points out that
echo # line 114 here: p *= zk + w * z * (zk - 1.0)/(z-1.0);
echo # is wrong, and therefore unfairly disparages Tolman's improvement.
echo # Gary Perlman <perlman@cis.ohio-state.edu> agrees, but says even
echo # the slow version is plenty fast enough, and since this version
echo # works correctly, he is inclined to leave well enough alone.
echo ########
# This is a shell archive, meaning:
# 1. Remove everything above the #! /bin/sh line.
# 2. Save the resulting text in a file.
# 3. Execute the file with /bin/sh (not csh) to create:
# z.c
# chisq.c
# f.c
# This archive created: Wed May 27 02:19:35 1987
# By: Gary Perlman (Wang Institute, Tyngsboro, MA 01879 USA)
export PATH; PATH=/bin:/usr/bin:$PATH
echo shar: "extracting 'z.c'" '(3123 characters)'
if test -f 'z.c'
then
echo shar: "will not over-write existing file 'z.c'"
else
sed 's/^ X//' << \SHAR_EOF > 'z.c'
X/*HEADER
X Module: z.c
X Purpose: compute approximations to normal z distribution probabilities
X Programmer: Gary Perlman
X Organization: Wang Institute, Tyngsboro, MA 01879
X Tester: compile with -DZTEST to include main program
X Copyright: none
X Tabstops: 4
X*/
X
X/*LINTLIBRARY*/
Xstatic char sccsfid[] = "@(#) z.c 5.1 (|stat) 12/26/85";
X#include <math.h>
X
X#define Z_EPSILON 0.000001 /* accuracy of critz approximation */
X#define Z_MAX 6.0 /* maximum meaningful z value */
X
Xdouble poz ();
Xdouble critz ();
X
X#ifdef ZTEST
Xmain ()
X {
X double z;
X printf ("%4s %10s %10s %10s\n",
X "z", "poz(z)", "poz(-z)", "z'");
X for (z = 0.0; z <= Z_MAX; z += .01)
X {
X printf ("%4.2f %10.6f %10.6f %10.6f\n",
X z, poz (z), poz (-z), critz (poz (z)));
X }
X }
X#endif ZTEST
X
X /*FUNCTION poz: probability of normal z value */
X/*ALGORITHM
X Adapted from a polynomial approximation in:
X Ibbetson D, Algorithm 209
X Collected Algorithms of the CACM 1963 p. 616
X Note:
X This routine has six digit accuracy, so it is only useful for absolute
X z values < 6. For z values >= to 6.0, poz() returns 0.0.
X*/
Xdouble /*VAR returns cumulative probability from -oo to z */
Xpoz (z)
Xdouble z; /*VAR normal z value */
X {
X double y, x, w;
X
X if (z == 0.0)
X x = 0.0;
X else
X {
X y = 0.5 * fabs (z);
X if (y >= (Z_MAX * 0.5))
X x = 1.0;
X else if (y < 1.0)
X {
X w = y*y;
X x = ((((((((0.000124818987 * w
X -0.001075204047) * w +0.005198775019) * w
X -0.019198292004) * w +0.059054035642) * w
X -0.151968751364) * w +0.319152932694) * w
X -0.531923007300) * w +0.797884560593) * y * 2.0;
X }
X else
X {
X y -= 2.0;
X x = (((((((((((((-0.000045255659 * y
X +0.000152529290) * y -0.000019538132) * y
X -0.000676904986) * y +0.001390604284) * y
X -0.000794620820) * y -0.002034254874) * y
X +0.006549791214) * y -0.010557625006) * y
X +0.011630447319) * y -0.009279453341) * y
X +0.005353579108) * y -0.002141268741) * y
X +0.000535310849) * y +0.999936657524;
X }
X }
X return (z > 0.0 ? ((x + 1.0) * 0.5) : ((1.0 - x) * 0.5));
X }
X
X /*FUNCTION critz: compute critical z value to produce given probability */
X/*ALGORITHM
X Begin with upper and lower limits for z values (maxz and minz)
X set to extremes. Choose a z value (zval) between the extremes.
X Compute the probability of the z value. Set minz or maxz, based
X on whether the probability is less than or greater than the
X desired p. Continue adjusting the extremes until they are
X within Z_EPSILON of each other.
X*/
Xdouble /*VAR returns z such that fabs (poz(p) - z) <= .000001 */
Xcritz (p)
Xdouble p; /*VAR critical probability level */
X {
X double minz = -Z_MAX; /* minimum of range of z */
X double maxz = Z_MAX; /* maximum of range of z */
X double zval = 0.0; /* computed/returned z value */
X double poz (), pval; /* prob (z) function, pval := poz (zval) */
X
X if (p <= 0.0 || p >= 1.0)
X return (0.0);
X
X while (maxz - minz > Z_EPSILON)
X {
X pval = poz (zval);
X if (pval > p)
X maxz = zval;
X else
X minz = zval;
X zval = (maxz + minz) * 0.5;
X }
X return (zval);
X }
SHAR_EOF
if test 3123 -ne "`wc -c < 'z.c'`"
then
echo shar: "error transmitting 'z.c'" '(should have been 3123 characters)'
fi
fi
echo shar: "extracting 'chisq.c'" '(3003 characters)'
if test -f 'chisq.c'
then
echo shar: "will not over-write existing file 'chisq.c'"
else
sed 's/^ X//' << \SHAR_EOF > 'chisq.c'
X/*
X Module: chisq.c
X Purpose: compute approximations to chisquare distribution probabilities
X Contents: pochisq(), critchi()
X Uses: poz() in z.c (Algorithm 209)
X Programmer: Gary Perlman
X Organization: Wang Institute, Tyngsboro, MA 01879
X Tester: compile with -DCHISQTEST to include main program
X Copyright: none
X Tabstops: 4
X*/
X
X/*LINTLIBRARY*/
X#include <math.h>
Xstatic char sccsfid[] = "@(#) chisq.c 5.2 (|stat) 12/08/86";
X
X#define CHI_EPSILON 0.000001 /* accuracy of critchi approximation */
X#define CHI_MAX 99999.0 /* maximum chi square value */
X
X#define LOG_SQRT_PI 0.5723649429247000870717135 /* log (sqrt (pi)) */
X#define I_SQRT_PI 0.5641895835477562869480795 /* 1 / sqrt (pi) */
X#define BIGX 20.0 /* max value to represent exp (x) */
X#define ex(x) (((x) < -BIGX) ? 0.0 : exp (x))
X
Xdouble pochisq ();
Xdouble critchi ();
X
X#ifdef CHISQTEST
Xdouble Prob[] = { .10, .05, .01, .005, .001, -1.0 };
Xmain ()
X {
X int df;
X int p;
X printf ("%-4s ", "df");
X for (p = 0; Prob[p] > 0.0; p++)
X printf ("%8.3f ", Prob[p]);
X putchar ('\n');
X for (df = 1; df < 30; df++)
X {
X printf ("%4d ", df);
X for (p = 0; Prob[p] > 0.0; p++)
X printf ("%8.3f ", critchi (Prob[p], df));
X putchar ('\n');
X }
X }
X#endif CHISQTEST
X
X /*FUNCTION pochisq: probability of chi sqaure value */
X/*ALGORITHM Compute probability of chi square value.
X Adapted from:
X Hill, I. D. and Pike, M. C. Algorithm 299
X Collected Algorithms for the CACM 1967 p. 243
X Updated for rounding errors based on remark in
X ACM TOMS June 1985, page 185
X*/
Xdouble
Xpochisq (x, df)
Xdouble x; /* obtained chi-square value */
Xint df; /* degrees of freedom */
X {
X double a, y, s;
X double e, c, z;
X double poz (); /* computes probability of normal z score */
X int even; /* true if df is an even number */
X
X if (x <= 0.0 || df < 1)
X return (1.0);
X
X a = 0.5 * x;
X even = (2*(df/2)) == df;
X if (df > 1)
X y = ex (-a);
X s = (even ? y : (2.0 * poz (-sqrt (x))));
X if (df > 2)
X {
X x = 0.5 * (df - 1.0);
X z = (even ? 1.0 : 0.5);
X if (a > BIGX)
X {
X e = (even ? 0.0 : LOG_SQRT_PI);
X c = log (a);
X while (z <= x)
X {
X e = log (z) + e;
X s += ex (c*z-a-e);
X z += 1.0;
X }
X return (s);
X }
X else
X {
X e = (even ? 1.0 : (I_SQRT_PI / sqrt (a)));
X c = 0.0;
X while (z <= x)
X {
X e = e * (a / z);
X c = c + e;
X z += 1.0;
X }
X return (c * y + s);
X }
X }
X else
X return (s);
X }
X
X /*FUNCTION critchi: compute critical chi square value to produce given p */
Xdouble
Xcritchi (p, df)
Xdouble p;
Xint df;
X {
X double minchisq = 0.0;
X double maxchisq = CHI_MAX;
X double chisqval;
X
X if (p <= 0.0)
X return (maxchisq);
X else if (p >= 1.0)
X return (0.0);
X
X chisqval = df / sqrt (p); /* fair first value */
X while (maxchisq - minchisq > CHI_EPSILON)
X {
X if (pochisq (chisqval, df) < p)
X maxchisq = chisqval;
X else
X minchisq = chisqval;
X chisqval = (maxchisq + minchisq) * 0.5;
X }
X return (chisqval);
X }
SHAR_EOF
if test 3003 -ne "`wc -c < 'chisq.c'`"
then
echo shar: "error transmitting 'chisq.c'" '(should have been 3003 characters)'
fi
fi
echo shar: "extracting 'f.c'" '(3847 characters)'
if test -f 'f.c'
then
echo shar: "will not over-write existing file 'f.c'"
else
sed 's/^ X//' << \SHAR_EOF > 'f.c'
X/*
X Module: f.c
X Purpose: compute approximations to F distribution probabilities
X Contents: pof(), critf()
X Programmer: Gary Perlman
X Organization: Wang Institute, Tyngsboro, MA 01879
X Tester: compile with -DFTEST to include main program
X Copyright: none
X Tabstops: 4
X*/
X
X/*LINTLIBRARY*/
X#include <math.h>
Xstatic char sccsfid[] = "@(#) f.c 5.2 (|stat) 12/26/85";
X
X#ifndef I_PI /* 1 / pi */
X#define I_PI 0.3183098861837906715377675
X#endif
X#define F_EPSILON 0.000001 /* accuracy of critf approximation */
X#define F_MAX 9999.0 /* maximum F ratio */
X
Xdouble pof ();
Xdouble critf ();
X
X#ifdef FTEST
X
Xint DFs[] = { 1, 2, 5, 10, 20, 40, 60, 120, -1 };
Xdouble Prob[] = { .10, .05, .01, .005, .001, -1.0 };
X
Xmain ()
X {
X int dfnum;
X int dfdenom;
X int p;
X
X for (p = 0; Prob[p] > 0.0; p++)
X {
X printf ("alpha = %g df1\n", Prob[p]);
X printf ("%-4s\\", "df2");
X for (dfnum = 0; DFs[dfnum] > 0; dfnum++)
X printf ("%7d ", DFs[dfnum]);
X putchar ('\n');
X for (dfdenom = 0; DFs[dfdenom] > 0; dfdenom++)
X {
X printf ("%4d ", DFs[dfdenom]);
X for (dfnum = 0; DFs[dfnum] > 0; dfnum++)
X printf ("%7.2f ", critf (Prob[p], DFs[dfnum], DFs[dfdenom]));
X putchar ('\n');
X }
X putchar ('\n');
X }
X }
X#endif FTEST
X
X /*FUNCTION pof: probability of F */
X/*ALGORITHM Compute probability of F ratio.
X Adapted from Collected Algorithms of the CACM
X Algorithm 322
X Egon Dorrer
X*/
Xdouble
Xpof (F, df1, df2)
Xdouble F;
Xint df1, df2;
X {
X int i, j;
X int a, b;
X double w, y, z, d, p;
X
X if (F < F_EPSILON || df1 < 1 || df2 < 1)
X return (1.0);
X a = df1%2 ? 1 : 2;
X b = df2%2 ? 1 : 2;
X w = (F * df1) / df2;
X z = 1.0 / (1.0 + w);
X if (a == 1)
X if (b == 1)
X {
X p = sqrt (w);
X y = I_PI; /* 1 / 3.14159 */
X d = y * z / p;
X p = 2.0 * y * atan (p);
X }
X else
X {
X p = sqrt (w * z);
X d = 0.5 * p * z / w;
X }
X else if (b == 1)
X {
X p = sqrt (z);
X d = 0.5 * z * p;
X p = 1.0 - p;
X }
X else
X {
X d = z * z;
X p = w * z;
X }
X y = 2.0 * w / z;
X#ifdef REMARK /* speedup modification suggested by Tolman (wrong answer!) */
X if (a == 1)
X for (j = b + 2; j <= df2; j += 2)
X {
X d *= (1.0 + a / (j - 2.0)) * z;
X p += d * y / (j - 1.0);
X }
X else
X {
X double zk = 1.0;
X for (j = (df2 - 1) / 2; j; j--)
X zk *= z;
X d *= zk * df2/b;
X p *= zk + w * z * (zk - 1.0)/(z-1.0);
X }
X#else /* original version */
X for (j = b + 2; j <= df2; j += 2)
X {
X d *= (1.0 + a / (j - 2.0)) * z;
X p = (a == 1 ? p + d * y / (j - 1.0) : (p + w) * z);
X }
X#endif REMARK
X y = w * z;
X z = 2.0 / z;
X b = df2 - 2;
X for (i = a + 2; i <= df1; i += 2)
X {
X j = i + b;
X d *= y * j / (i - 2.0);
X p -= z * d / j;
X }
X /* correction for approximation errors suggested in certification */
X if (p < 0.0)
X p = 0.0;
X else if (p > 1.0)
X p = 1.0;
X return (1.0-p);
X }
X
X /*FUNCTION critf: compute critical F value t produce given probability */
X/*ALGORITHM
X Begin with upper and lower limits for F values (maxf and minf)
X set to extremes. Choose an f value (fval) between the extremes.
X Compute the probability of the f value. Set minf or maxf, based
X on whether the probability is less than or greater than the
X desired p. Continue adjusting the extremes until they are
X within F_EPSILON of each other.
X*/
Xdouble
Xcritf (p, df1, df2)
Xdouble p;
Xint df1;
Xint df2;
X {
X double fval;
X double fabs (); /* floating absolute value */
X double maxf = F_MAX; /* maximum F ratio */
X double minf = 0.0; /* minimum F ratio */
X
X if (p <= 0.0 || p >= 1.0)
X return (0.0);
X
X fval = 1.0 / p; /* the smaller the p, the larger the F */
X
X while (fabs (maxf - minf) > F_EPSILON)
X {
X if (pof (fval, df1, df2) < p) /* F too large */
X maxf = fval;
X else /* F too small */
X minf = fval;
X fval = (maxf + minf) * 0.5;
X }
X
X return (fval);
X }
SHAR_EOF
if test 3847 -ne "`wc -c < 'f.c'`"
then
echo shar: "error transmitting 'f.c'" '(should have been 3847 characters)'
fi
fi
exit 0
# End of shell archive
moved to a/perlman