6 extern double gamma(), ppnd();
7 static double gammp(), gser(), gcf(), gnorm(), ppchi2();
15 double ppgamma(p
, a
, ifault
)
23 x
= ppchi2(&p
, &v
, &g
, ifault
);
32 From Numerical Recipes, with normal approximation from Appl. Stat. 239
35 #define EPSILON 1.0e-14
36 #define LARGE_A 10000.0
39 static double gammp(a
, x
)
44 if (x
<= 0.0 || a
<= 0.0) p
= 0.0;
45 else if (a
> LARGE_A
) p
= gnorm(a
, x
);
48 if (x
< (a
+ 1.0)) p
= gser(a
, x
, gln
);
49 else p
= 1.0 - gcf(a
, x
, gln
);
54 /* compute gamma cdf by a normal approximation */
55 static double gnorm(a
, x
)
60 if (x
<= 0.0 || a
<= 0.0) p
= 0.0;
62 sx
= sqrt(a
) * 3.0 * (pow(x
/ a
, 1.0 / 3.0) + 1.0 / (a
* 9.0) - 1.0);
68 /* compute gamma cdf by its series representation */
69 static double gser(a
, x
, gln
)
72 double p
, sum
, del
, ap
;
75 if (x
<= 0.0 || a
<= 0.0) p
= 0.0;
79 for (n
= 1; ! done
&& n
< ITMAX
; n
++) {
83 if (fabs(del
) < EPSILON
) done
= TRUE
;
85 p
= sum
* exp(- x
+ a
* log(x
) - gln
);
90 /* compute complementary gamma cdf by its continued fraction expansion */
91 static double gcf(a
, x
, gln
)
94 double gold
= 0.0, g
, fac
= 1.0, b1
= 1.0;
95 double b0
= 0.0, anf
, ana
, an
, a1
, a0
= 1.0;
101 for(an
= 1.0; ! done
&& an
<= ITMAX
; an
+= 1.0) {
103 a0
= (a1
+ a0
* ana
) * fac
;
104 b0
= (b1
+ b0
* ana
) * fac
;
106 a1
= x
* a0
+ anf
* a1
;
107 b1
= x
* b0
+ anf
* b1
;
111 if (fabs((g
- gold
) / g
) < EPSILON
) {
112 p
= exp(-x
+ a
* log(x
) - gln
) * g
;
121 static double gammad(x
, a
, iflag
)
127 gammabase(x
, a
, &cdf
);
132 ppchi2.f -- translated by f2c and modified
134 Algorithm AS 91 Appl. Statist. (1975) Vol.24, P.35
135 To evaluate the percentage points of the chi-squared
136 probability distribution function.
138 p must lie in the range 0.000002 to 0.999998,
139 (but I am using it for 0 < p < 1 - seems to work)
141 g must be supplied and should be equal to ln(gamma(v/2.0))
143 Auxiliary routines required: ppnd = AS 111 (or AS 241) and gammad.
146 static double ppchi2(p
, v
, g
, ifault
)
150 /* Initialized data */
152 static double aa
= .6931471806;
153 static double six
= 6.;
154 static double c1
= .01;
155 static double c2
= .222222;
156 static double c3
= .32;
157 static double c4
= .4;
158 static double c5
= 1.24;
159 static double c6
= 2.2;
160 static double c7
= 4.67;
161 static double c8
= 6.66;
162 static double c9
= 6.73;
163 static double e
= 5e-7;
164 static double c10
= 13.32;
165 static double c11
= 60.;
166 static double c12
= 70.;
167 static double c13
= 84.;
168 static double c14
= 105.;
169 static double c15
= 120.;
170 static double c16
= 127.;
171 static double c17
= 140.;
172 static double c18
= 1175.;
173 static double c19
= 210.;
174 static double c20
= 252.;
175 static double c21
= 2264.;
176 static double c22
= 294.;
177 static double c23
= 346.;
178 static double c24
= 420.;
179 static double c25
= 462.;
180 static double c26
= 606.;
181 static double c27
= 672.;
182 static double c28
= 707.;
183 static double c29
= 735.;
184 static double c30
= 889.;
185 static double c31
= 932.;
186 static double c32
= 966.;
187 static double c33
= 1141.;
188 static double c34
= 1182.;
189 static double c35
= 1278.;
190 static double c36
= 1740.;
191 static double c37
= 2520.;
192 static double c38
= 5040.;
193 static double zero
= 0.;
194 static double half
= .5;
195 static double one
= 1.;
196 static double two
= 2.;
197 static double three
= 3.;
200 static double pmin = 2e-6;
201 static double pmax = .999998;
203 static double pmin
= 0.0;
204 static double pmax
= 1.0;
206 /* System generated locals */
207 double ret_val
, d_1
, d_2
;
209 /* Local variables */
210 static double a
, b
, c
, q
, t
, x
, p1
, p2
, s1
, s2
, s3
, s4
, s5
, s6
, ch
;
215 /* test arguments and initialise */
218 if (*p
<= pmin
|| *p
>= pmax
) return ret_val
;
220 if (*v
<= zero
) return ret_val
;
225 if (*v
< -c5
* log(*p
)) {
226 /* starting approximation for small chi-squared */
227 ch
= pow(*p
* xx
* exp(*g
+ xx
* aa
), one
/ xx
);
234 /* call to algorithm AS 111 - note that p has been tested above. */
235 /* AS 241 could be used as an alternative. */
238 /* starting approximation using Wilson and Hilferty estimate */
240 /* Computing 3rd power */
241 d_1
= x
* sqrt(p1
) + one
- p1
;
242 ch
= *v
* (d_1
* d_1
* d_1
);
244 /* starting approximation for p tending to 1 */
245 if (ch
> c6
* *v
+ six
)
246 ch
= -two
* (log(one
- *p
) - c
* log(half
* ch
) + *g
);
249 /* starting approximation for v less than or equal to 0.32 */
254 p1
= one
+ ch
* (c7
+ ch
);
255 p2
= ch
* (c9
+ ch
* (c8
+ ch
));
256 d_1
= -half
+ (c7
+ two
* ch
) / p1
;
257 d_2
= (c9
+ ch
* (c10
+ three
* ch
)) / p2
;
259 ch
-= (one
- exp(a
+ *g
+ half
* ch
+ c
* aa
) * p2
/ p1
) / t
;
260 } while (fabs(q
/ ch
- one
) > c1
);
264 /* call to gammad and calculation of seven term Taylor series */
267 p2
= *p
- gammad(&p1
, &xx
, &if1
);
272 t
= p2
* exp(xx
* aa
+ *g
+ p1
- c
* log(ch
));
274 a
= half
* t
- b
* c
;
275 s1
= (c19
+ a
* (c17
+ a
* (c14
+ a
* (c13
+ a
* (c12
+ c11
* a
))))) / c24
;
276 s2
= (c24
+ a
* (c29
+ a
* (c32
+ a
* (c33
+ c35
* a
)))) / c37
;
277 s3
= (c19
+ a
* (c25
+ a
* (c28
+ c31
* a
))) / c37
;
278 s4
= (c20
+ a
* (c27
+ c34
* a
) + c
* (c22
+ a
* (c30
+ c36
* a
))) / c38
;
279 s5
= (c13
+ c21
* a
+ c
* (c18
+ c26
* a
)) / c37
;
280 s6
= (c15
+ c
* (c23
+ c16
* c
)) / c38
;
281 d_1
= (s3
- b
* (s4
- b
* (s5
- b
* s6
)));
282 d_1
= (s1
- b
* (s2
- b
* d_1
));
283 ch
+= t
* (one
+ half
* t
* s1
- b
* c
* d_1
);
284 } while (fabs(q
/ ch
- one
) > e
);