6 extern void normbase(double *,double *);
7 extern double gamma(), ppnd();
8 static double gammp(), gser(), gcf(), gnorm(), ppchi2();
12 gammabase(double *x
, double *a
, double *p
)
18 ppgamma(double p
, double a
, int *ifault
)
24 x
= ppchi2(&p
, &v
, &g
, ifault
);
33 From Numerical Recipes, with normal approximation from Appl. Stat. 239
36 #define EPSILON 1.0e-14
37 #define LARGE_A 10000.0
40 static double gammp(a
, x
)
45 if (x
<= 0.0 || a
<= 0.0) p
= 0.0;
46 else if (a
> LARGE_A
) p
= gnorm(a
, x
);
49 if (x
< (a
+ 1.0)) p
= gser(a
, x
, gln
);
50 else p
= 1.0 - gcf(a
, x
, gln
);
55 /* compute gamma cdf by a normal approximation */
56 static double gnorm(a
, x
)
61 if (x
<= 0.0 || a
<= 0.0) p
= 0.0;
63 sx
= sqrt(a
) * 3.0 * (pow(x
/ a
, 1.0 / 3.0) + 1.0 / (a
* 9.0) - 1.0);
69 /* compute gamma cdf by its series representation */
70 static double gser(a
, x
, gln
)
73 double p
, sum
, del
, ap
;
76 if (x
<= 0.0 || a
<= 0.0) p
= 0.0;
80 for (n
= 1; ! done
&& n
< ITMAX
; n
++) {
84 if (fabs(del
) < EPSILON
) done
= TRUE
;
86 p
= sum
* exp(- x
+ a
* log(x
) - gln
);
91 /* compute complementary gamma cdf by its continued fraction expansion */
92 static double gcf(a
, x
, gln
)
95 double gold
= 0.0, g
, fac
= 1.0, b1
= 1.0;
96 double b0
= 0.0, anf
, ana
, an
, a1
, a0
= 1.0;
102 for(an
= 1.0; ! done
&& an
<= ITMAX
; an
+= 1.0) {
104 a0
= (a1
+ a0
* ana
) * fac
;
105 b0
= (b1
+ b0
* ana
) * fac
;
107 a1
= x
* a0
+ anf
* a1
;
108 b1
= x
* b0
+ anf
* b1
;
112 if (fabs((g
- gold
) / g
) < EPSILON
) {
113 p
= exp(-x
+ a
* log(x
) - gln
) * g
;
122 static double gammad(x
, a
, iflag
)
128 gammabase(x
, a
, &cdf
);
133 ppchi2.f -- translated by f2c and modified
135 Algorithm AS 91 Appl. Statist. (1975) Vol.24, P.35
136 To evaluate the percentage points of the chi-squared
137 probability distribution function.
139 p must lie in the range 0.000002 to 0.999998,
140 (but I am using it for 0 < p < 1 - seems to work)
142 g must be supplied and should be equal to ln(gamma(v/2.0))
144 Auxiliary routines required: ppnd = AS 111 (or AS 241) and gammad.
147 static double ppchi2(p
, v
, g
, ifault
)
151 /* Initialized data */
153 static double aa
= .6931471806;
154 static double six
= 6.;
155 static double c1
= .01;
156 static double c2
= .222222;
157 static double c3
= .32;
158 static double c4
= .4;
159 static double c5
= 1.24;
160 static double c6
= 2.2;
161 static double c7
= 4.67;
162 static double c8
= 6.66;
163 static double c9
= 6.73;
164 static double e
= 5e-7;
165 static double c10
= 13.32;
166 static double c11
= 60.;
167 static double c12
= 70.;
168 static double c13
= 84.;
169 static double c14
= 105.;
170 static double c15
= 120.;
171 static double c16
= 127.;
172 static double c17
= 140.;
173 static double c18
= 1175.;
174 static double c19
= 210.;
175 static double c20
= 252.;
176 static double c21
= 2264.;
177 static double c22
= 294.;
178 static double c23
= 346.;
179 static double c24
= 420.;
180 static double c25
= 462.;
181 static double c26
= 606.;
182 static double c27
= 672.;
183 static double c28
= 707.;
184 static double c29
= 735.;
185 static double c30
= 889.;
186 static double c31
= 932.;
187 static double c32
= 966.;
188 static double c33
= 1141.;
189 static double c34
= 1182.;
190 static double c35
= 1278.;
191 static double c36
= 1740.;
192 static double c37
= 2520.;
193 static double c38
= 5040.;
194 static double zero
= 0.;
195 static double half
= .5;
196 static double one
= 1.;
197 static double two
= 2.;
198 static double three
= 3.;
201 static double pmin = 2e-6;
202 static double pmax = .999998;
204 static double pmin
= 0.0;
205 static double pmax
= 1.0;
207 /* System generated locals */
208 double ret_val
, d_1
, d_2
;
210 /* Local variables */
211 static double a
, b
, c
, q
, t
, x
, p1
, p2
, s1
, s2
, s3
, s4
, s5
, s6
, ch
;
216 /* test arguments and initialise */
219 if (*p
<= pmin
|| *p
>= pmax
) return ret_val
;
221 if (*v
<= zero
) return ret_val
;
226 if (*v
< -c5
* log(*p
)) {
227 /* starting approximation for small chi-squared */
228 ch
= pow(*p
* xx
* exp(*g
+ xx
* aa
), one
/ xx
);
235 /* call to algorithm AS 111 - note that p has been tested above. */
236 /* AS 241 could be used as an alternative. */
239 /* starting approximation using Wilson and Hilferty estimate */
241 /* Computing 3rd power */
242 d_1
= x
* sqrt(p1
) + one
- p1
;
243 ch
= *v
* (d_1
* d_1
* d_1
);
245 /* starting approximation for p tending to 1 */
246 if (ch
> c6
* *v
+ six
)
247 ch
= -two
* (log(one
- *p
) - c
* log(half
* ch
) + *g
);
250 /* starting approximation for v less than or equal to 0.32 */
255 p1
= one
+ ch
* (c7
+ ch
);
256 p2
= ch
* (c9
+ ch
* (c8
+ ch
));
257 d_1
= -half
+ (c7
+ two
* ch
) / p1
;
258 d_2
= (c9
+ ch
* (c10
+ three
* ch
)) / p2
;
260 ch
-= (one
- exp(a
+ *g
+ half
* ch
+ c
* aa
) * p2
/ p1
) / t
;
261 } while (fabs(q
/ ch
- one
) > c1
);
265 /* call to gammad and calculation of seven term Taylor series */
268 p2
= *p
- gammad(&p1
, &xx
, &if1
);
273 t
= p2
* exp(xx
* aa
+ *g
+ p1
- c
* log(ch
));
275 a
= half
* t
- b
* c
;
276 s1
= (c19
+ a
* (c17
+ a
* (c14
+ a
* (c13
+ a
* (c12
+ c11
* a
))))) / c24
;
277 s2
= (c24
+ a
* (c29
+ a
* (c32
+ a
* (c33
+ c35
* a
)))) / c37
;
278 s3
= (c19
+ a
* (c25
+ a
* (c28
+ c31
* a
))) / c37
;
279 s4
= (c20
+ a
* (c27
+ c34
* a
) + c
* (c22
+ a
* (c30
+ c36
* a
))) / c38
;
280 s5
= (c13
+ c21
* a
+ c
* (c18
+ c26
* a
)) / c37
;
281 s6
= (c15
+ c
* (c23
+ c16
* c
)) / c38
;
282 d_1
= (s3
- b
* (s4
- b
* (s5
- b
* s6
)));
283 d_1
= (s1
- b
* (s2
- b
* d_1
));
284 ch
+= t
* (one
+ half
* t
* s1
- b
* c
* d_1
);
285 } while (fabs(q
/ ch
- one
) > e
);