DataMuseum.dk

top - download

    Length: 1865 (0x749)
    Types: TextFile
    Names: »chisq.c«

└─⟦a0efdde77⟧ Bits:30001252 EUUGD11 Tape, 1987 Spring Conference Helsinki
    └─ ⟦this⟧ »EUUGD11/stat-5.3/eu/stat/src/chisq.c« 

TextFile

/*  Copyright 1986 Gary Perlman */

/*LINTLIBRARY*/
#include "stat.h"
FUN(chisq,chi-square distribution functions,5.2,12/08/86)

#define	LOG_SQRT_PI       0.5723649429247000870717135 /* log (sqrt (pi)) */
#define	I_SQRT_PI         0.5641895835477562869480795 /* 1 / sqrt (pi) */
#define	BIGX              20.0   /* max value to represent exp (x) */
#define	ex(x)             (((x) < -BIGX) ? 0.0 : exp (x))

/*
	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
*/
double
pochisq (x, df)
double	x;      /* obtained chi-square value */
int 	df;      /* degrees of freedom */
	{
	double	a, y, s;
	double	poz ();   /* computes probability of normal z score */
	int 	even;

	if (x <= 0.0 || df < 1)
		return (1.0);

	a = 0.5 * x;
	even = (2*(df/2)) == df;
	if (df > 1)
		y = ex (-a);
	s = (even ? y : (2.0 * poz (-sqrt (x))));
	if (df > 2)
		{
		double	e, c, z;
		x = 0.5 * (df - 1.0);
		z = (even ? 1.0 : 0.5);
		if (a > BIGX)
			{
			e = (even ? 0.0 : LOG_SQRT_PI);
			c = log (a);
			while (z <= x)
				{
				e = log (z) + e;
				s += ex (c*z-a-e);
				z += 1.0;
				}
			return (s);
			}
		else
			{
			e = (even ? 1.0 : (I_SQRT_PI / sqrt (a)));
			c = 0.0;
			while (z <= x)
				{
				e = e * (a / z);
				c = c + e;
				z += 1.0;
				}
			return (c * y + s);
			}
		}
	else
		return (s);
	}

double
critchi (p, df)
double	p;
int 	df;
	{
	double	minchisq = 0.0;
	double	maxchisq = 99999.0;
	double	chisqval;
	double	pochisq ();

	if (p <= 0.0)
		return (maxchisq);
	else if (p >= 1.0)
		return (0.0);

	chisqval = 1.0 / p;
	while (maxchisq - minchisq > .000001)
		{
		if (pochisq (chisqval, df) < p)
			maxchisq = chisqval;
		else
			minchisq = chisqval;
		chisqval = (maxchisq + minchisq) * 0.5;
		}
	return (chisqval);
	}