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
137
138
139
140
141
142
143
144
|
/* This code was taken from http://www.fourmilab.ch/random/ where it states that:
This software is in the public domain. Permission to use, copy, modify, and distribute
this software and its documentation for any purpose and without fee is hereby granted,
without any conditions or restrictions. This software is provided “as is” without
express or implied warranty. */
/*
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;
y = 0.0;
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;
}
}
|
