clean up of consistent-dataframe-p generic/methods. They were subtly wrong.
[CommonLispStat.git] / lib / gammabase.c
blob66f833e5abd45d6646d476c219346c35853a91ef
1 #include "xmath.h"
3 #define TRUE 1
4 #define FALSE 0
6 extern void normbase(double *,double *);
7 extern double gamma(), ppnd();
8 static double gammp(), gser(), gcf(), gnorm(), ppchi2();
11 void
12 gammabase(double *x, double *a, double *p)
14 *p = gammp(*a, *x);
17 double
18 ppgamma(double p, double a, int *ifault)
20 double x, v, g;
22 v = 2.0 * a;
23 g = gamma(a);
24 x = ppchi2(&p, &v, &g, ifault);
25 return(x / 2.0);
29 Static Routines
33 From Numerical Recipes, with normal approximation from Appl. Stat. 239
36 #define EPSILON 1.0e-14
37 #define LARGE_A 10000.0
38 #define ITMAX 1000
40 static double gammp(a, x)
41 double a, x;
43 double gln, p;
45 if (x <= 0.0 || a <= 0.0) p = 0.0;
46 else if (a > LARGE_A) p = gnorm(a, x);
47 else {
48 gln = gamma(a);
49 if (x < (a + 1.0)) p = gser(a, x, gln);
50 else p = 1.0 - gcf(a, x, gln);
52 return(p);
55 /* compute gamma cdf by a normal approximation */
56 static double gnorm(a, x)
57 double a, x;
59 double p, sx;
61 if (x <= 0.0 || a <= 0.0) p = 0.0;
62 else {
63 sx = sqrt(a) * 3.0 * (pow(x / a, 1.0 / 3.0) + 1.0 / (a * 9.0) - 1.0);
64 normbase(&sx, &p);
66 return(p);
69 /* compute gamma cdf by its series representation */
70 static double gser(a, x, gln)
71 double a, x, gln;
73 double p, sum, del, ap;
74 int n, done = FALSE;
76 if (x <= 0.0 || a <= 0.0) p = 0.0;
77 else {
78 ap = a;
79 del = sum = 1.0 / a;
80 for (n = 1; ! done && n < ITMAX; n++) {
81 ap += 1.0;
82 del *= x / ap;
83 sum += del;
84 if (fabs(del) < EPSILON) done = TRUE;
86 p = sum * exp(- x + a * log(x) - gln);
88 return(p);
91 /* compute complementary gamma cdf by its continued fraction expansion */
92 static double gcf(a, x, gln)
93 double 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;
97 double p;
98 int done = FALSE;
100 a1 = x;
101 p = 0.0;
102 for(an = 1.0; ! done && an <= ITMAX; an += 1.0) {
103 ana = an - a;
104 a0 = (a1 + a0 * ana) * fac;
105 b0 = (b1 + b0 * ana) * fac;
106 anf = an * fac;
107 a1 = x * a0 + anf * a1;
108 b1 = x * b0 + anf * b1;
109 if (a1 != 0.0) {
110 fac = 1.0 / a1;
111 g = b1 * fac;
112 if (fabs((g - gold) / g) < EPSILON) {
113 p = exp(-x + a * log(x) - gln) * g;
114 done = TRUE;
116 gold = g;
119 return(p);
122 static double gammad(x, a, iflag)
123 double *x, *a;
124 int *iflag;
126 double cdf;
128 gammabase(x, a, &cdf);
129 return(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)
141 v must be positive,
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)
148 double *p, *v, *g;
149 int *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;
212 static double xx;
213 static int if1;
216 /* test arguments and initialise */
217 ret_val = -one;
218 *ifault = 1;
219 if (*p <= pmin || *p >= pmax) return ret_val;
220 *ifault = 2;
221 if (*v <= zero) return ret_val;
222 *ifault = 0;
223 xx = half * *v;
224 c = xx - one;
226 if (*v < -c5 * log(*p)) {
227 /* starting approximation for small chi-squared */
228 ch = pow(*p * xx * exp(*g + xx * aa), one / xx);
229 if (ch < e) {
230 ret_val = ch;
231 return ret_val;
234 else if (*v > c3) {
235 /* call to algorithm AS 111 - note that p has been tested above. */
236 /* AS 241 could be used as an alternative. */
237 x = ppnd(*p, &if1);
239 /* starting approximation using Wilson and Hilferty estimate */
240 p1 = c2 / *v;
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);
249 else{
250 /* starting approximation for v less than or equal to 0.32 */
251 ch = c4;
252 a = log(one - *p);
253 do {
254 q = ch;
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;
259 t = d_1 - d_2;
260 ch -= (one - exp(a + *g + half * ch + c * aa) * p2 / p1) / t;
261 } while (fabs(q / ch - one) > c1);
264 do {
265 /* call to gammad and calculation of seven term Taylor series */
266 q = ch;
267 p1 = half * ch;
268 p2 = *p - gammad(&p1, &xx, &if1);
269 if (if1 != 0) {
270 *ifault = 3;
271 return ret_val;
273 t = p2 * exp(xx * aa + *g + p1 - c * log(ch));
274 b = t / 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);
287 ret_val = ch;
288 return ret_val;