blob: baee4131fe16decf01c41f0079b4bb10046b886f (
plain)
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
|
/*
ChiSquare.c
the chi-square cdf
after W.J. Kennedy and J.E. Gentle,
Statistical computing, p. 116
last modified 9 Feb 05 th
*/
#ifdef HAVE_ERF
#define Erf erf
#else
#include "Erf.c"
#endif
static inline real Normal(creal x)
{
return .5*Erf(x/1.414213562373095048801689) + .5;
}
/*********************************************************************/
static real ChiSquare(creal x, cint df)
{
real y;
if( df <= 0 ) return -999;
if( x <= 0 ) return 0;
if( x > 1000*df ) return 1;
if( df > 1000 ) {
if( x < 2 ) return 0;
y = 2./(9*df);
y = (pow(x/df, 1/3.) - (1 - y))/sqrt(y);
if( y > 5 ) return 1;
if( y < -18.8055 ) return 0;
return Normal(y);
}
y = .5*x;
if( df & 1 ) {
creal sqrty = sqrt(y);
real h = Erf(sqrty);
count i;
if( df == 1 ) return h;
y = sqrty*exp(-y)/.8862269254527579825931;
for( i = 3; i < df; i += 2 ) {
h -= y;
y *= x/i;
}
y = h - y;
}
else {
real term = exp(-y), sum = term;
count i;
for( i = 1; i < df/2; ++i )
sum += term *= y/i;
y = 1 - sum;
}
return Max(0., y);
}
|
