File: chisq.c

package info (click to toggle)
ent 1.2debian-2
  • links: PTS, VCS
  • area: main
  • in suites: bullseye, buster, sid
  • size: 260 kB
  • sloc: ansic: 410; makefile: 41
file content (136 lines) | stat: -rw-r--r-- 3,663 bytes parent folder | download | duplicates (2)
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
/*

    Compute probability of measured Chi Square value.
    
    This code was developed by Gary Perlman of the Wang
    Institute (full citation below) and has been minimally
    modified for use in this program.
    
*/

#include <math.h>

/*HEADER
	Module:       z.c
	Purpose:      compute approximations to normal z distribution probabilities
	Programmer:   Gary Perlman
	Organization: Wang Institute, Tyngsboro, MA 01879
	Copyright:    none
	Tabstops:     4
*/

#define	Z_MAX          6.0            /* maximum meaningful z value */

/*FUNCTION poz: probability of normal z value */
/*ALGORITHM
	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 >= to 6.0, poz() returns 0.0.
*/
static double        /*VAR returns cumulative probability from -oo to z */
poz(const double z)  /*VAR normal z value */
{
    double y, x, w;

    if (z == 0.0) {
    	x = 0.0;
    } else {
	y = 0.5 * fabs(z);
	if (y >= (Z_MAX * 0.5)) {
    	    x = 1.0;
	} else if (y < 1.0) {
	   w = y * y;
	   x = ((((((((0.000124818987 * w
		   -0.001075204047) * w +0.005198775019) * w
		   -0.019198292004) * w +0.059054035642) * w
		   -0.151968751364) * w +0.319152932694) * w
		   -0.531923007300) * w +0.797884560593) * y * 2.0;
	} else {
	    y -= 2.0;
	    x = (((((((((((((-0.000045255659 * y
		    +0.000152529290) * y -0.000019538132) * y
		    -0.000676904986) * y +0.001390604284) * y
		    -0.000794620820) * y -0.002034254874) * y
		    +0.006549791214) * y -0.010557625006) * y
		    +0.011630447319) * y -0.009279453341) * y
		    +0.005353579108) * y -0.002141268741) * y
		    +0.000535310849) * y +0.999936657524;
    	}
    }
    return (z > 0.0 ? ((x + 1.0) * 0.5) : ((1.0 - x) * 0.5));
}

/*
	Module:       chisq.c
	Purpose:      compute approximations to chisquare distribution probabilities
	Contents:     pochisq()
	Uses:         poz() in z.c (Algorithm 209)
	Programmer:   Gary Perlman
	Organization: Wang Institute, Tyngsboro, MA 01879
	Copyright:    none
	Tabstops:     4
*/

#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))

/*FUNCTION pochisq: probability of chi sqaure value */
/*ALGORITHM Compute probability of chi square value.
	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(
    	const double ax,    /* obtained chi-square value */
     	const int df	    /* degrees of freedom */
     	)
{
    double x = ax;
    double a, y, s;
    double e, c, z;
    int even;	    	    /* true if df is an even number */

    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) {
	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;
    }
}