1 /*
2 * $Id: gammp.i,v 1.1 2005-09-18 22:06:14 dhmunro Exp $
3 * Incomplete gamma (chi-square distribution) and beta functions.
4 */
5 /* Copyright (c) 2005, The Regents of the University of California.
6 * All rights reserved.
7 * This file is part of yorick (http://yorick.sourceforge.net).
8 * Read the accompanying LICENSE file for details.
9 */
10
11 /* Algorithms from Numerical Recipes, Press, et. al., section 6.2,
12 * Cambridge University Press (1988).
13 */
14
15 require, "gamma.i";
16
17 func gammp(a, x, &q, &lg)
18 /* DOCUMENT gammp(a, x)
19 * or gammp(a, x, q, lg)
20 * return P(a,x) = int[0 to x]{ du * u^(a-1)*exp(-u) } / gamma(a)
21 * optionally return Q(a,x) = 1-P(a,x) and ln(gamma(a))
22 *
23 * Note that erf(x)=gammp(0.5,x^2) and erfc(x)=gammq(0.5,x^2)
24 * Also P(chi2|nu)=gammp(0.5*nu,0.5*chi2)
25 * and Q(chi2|nu)=gammq(0.5*nu,0.5*chi2)
26 * are the probabilities that an observed chi-square be less than
27 * or greater than (P or Q) chi2 when there are nu degrees of freedom
28 * (terms in the chi-square sum).
29 *
30 * SEE ALSO: gammq, betai, ln_gamma (gamma.i), (dawson.i)
31 */
32 {
33 lg = ln_gamma(a);
34 fact = double(!x);
35 fact = (1.-fact)*exp(-x+a*log(x+fact)-lg);
36 x += 0.*a;
37 a += 0.*x;
38 mask = x < a+1.;
39 list = where(mask);
40 if (numberof(list)) {
41 p1 = fact(list)*_gammp(a(list), x(list));
42 q1 = 1. - p1;
43 }
44 list = where(!mask);
45 if (numberof(list)) {
46 q = fact(list)*_gammq(a(list), x(list));
47 p = 1. - q;
48 }
49 q = merge(q1, q, mask);
50 return merge(p1, p, mask);
51 }
52
53 func gammq(a, x, &p, &lg)
54 /* DOCUMENT gammq(a, x)
55 * or gammq(a, x, p, lg)
56 * return Q(a,x) = 1 - int[0 to x]{ du * u^(a-1)*exp(-u) } / gamma(a)
57 * optionally return P(a,x) = 1-Q(a,x) and ln(gamma(a))
58 *
59 * Note that erf(x)=gammp(0.5,x^2) and erfc(x)=gammq(0.5,x^2)
60 * Also P(chi2|nu)=gammp(0.5*nu,0.5*chi2)
61 * and Q(chi2|nu)=gammq(0.5*nu,0.5*chi2)
62 * are the probabilities that an observed chi-square be less than
63 * or greater than (P or Q) chi2 when there are nu degrees of freedom
64 * (terms in the chi-square sum).
65 *
66 * SEE ALSO: gammp, betai, ln_gamma (gamma.i), (dawson.i)
67 */
68 {
69 { local q; }
70 p = gammp(a, x, q, lg);
71 return q;
72 }
73
_gammp(a,x)74 func _gammp(a, x)
75 {
76 term = tot = total = 1./a;
77 ndx = indgen(numberof(total));
78 if (numberof(total) == 1) total = [total(1)];
79 for (i=0 ; i<_gamm_max ; ++i) {
80 a += 1.0;
81 term *= x/a;
82 tot += term;
83 mask = abs(term) < abs(tot)*_gamm_err;
84 list = where(mask);
85 if (numberof(list)) total(ndx(list)) = tot(list);
86 list = where(!mask);
87 if (!numberof(list)) break;
88 ndx = ndx(list);
89 a = a(list);
90 x = x(list);
91 term = term(list);
92 tot = tot(list);
93 }
94 if (numberof(list))
95 error, "global variable _gamm_err or _gamm_max too small";
96 if (numberof(total) == 1) total = total(1);
97 return total;
98 }
99
_gammq(a,x)100 func _gammq(a, x)
101 {
102 r = rold = b0 = 0.*x;
103 if (numberof(r) == 1) r = [r(1)];
104 a0 = b1 = ren = 1. + b0;
105 a1 = x;
106 ndx = indgen(numberof(r));
107 for (i=1 ; i<=_gamm_max ; ++i) {
108 tmp = i - a;
109 a0 = (a1+a0*tmp)*ren;
110 b0 = (b1+b0*tmp)*ren;
111 tmp = i*ren;
112 a1 = x*a0+tmp*a1;
113 b1 = x*b0+tmp*b1;
114 mask = abs(b1-a1*rold) < abs(b1)*_gamm_err;
115 list = where(mask);
116 if (numberof(list)) r(ndx(list)) = b1(list)/a1(list);
117 list = where(!mask);
118 if (!numberof(list)) break;
119 ndx = ndx(list);
120 a = a(list);
121 x = x(list);
122 a0 = a0(list);
123 a1 = a1(list);
124 b0 = b0(list);
125 b1 = b1(list);
126 ren = 1./(a1+!a1);
127 rold = b1*ren;
128 }
129 if (numberof(list))
130 error, "global variable _gamm_err or _gamm_max too small";
131 if (numberof(r) == 1) r = r(1);
132 return r;
133 }
134
135 _gamm_err = 1.e-13;
136 _gamm_max = 200;
137
betai(a,b,x)138 func betai(a, b, x)
139 /* DOCUMENT betai(a, b, x)
140 * return I_x(a,x) = int[0 to x]{ du * u^(a-1)*(1-u)^(b-1) } / beta(a,b)
141 * the incomplete beta function
142 *
143 * betai(a,b,x) = 1 - betai(b,a,1-x)
144 *
145 * Note that Student's t-distribution is
146 * A(t|nu) = 1 - betai(0.5*nu,0.5, nu/(nu+t^2))
147 * The F-distribution is
148 * Q(F|nu1,nu2) = betai(0.5*nu2,0.5*nu1, nu2/(nu2+F*nu1))
149 *
150 * SEE ALSO: gammp, gammq, ln_gamma (gamma.i), (dawson.i)
151 */
152 {
153 mask = x < (a+1.)/(a+b+2.);
154 fact = 0.*mask;
155 x += fact;
156 a += fact;
157 b += fact;
158 xc = 1.-x;
159 fact = double((!x) | (!xc));
160 fact = (1.-fact)*exp(ln_gamma(a+b)-ln_gamma(a)-ln_gamma(b)+
161 a*log(x+fact)+b*log(xc+fact));
162 list = where(mask);
163 if (numberof(list)) {
164 r1 = a(list);
165 r1 = fact(list)*_betai(r1, b(list), x(list))/r1;
166 }
167 list = where(!mask);
168 if (numberof(list)) {
169 r = b(list);
170 r = 1. - fact(list)*_betai(r, a(list), xc(list))/r;
171 }
172 return merge(r1, r, mask);
173 }
174
_betai(a,b,x)175 func _betai(a, b, x)
176 {
177 r = 0.*x;
178 a0 = b0 = b1 = rold = 1. + r;
179 a1 = 1. - (a+b)*x/(a+1.);
180 if (numberof(r) == 1) r = [r(1)];
181 ndx = indgen(numberof(r));
182 for (i=1 ; i<=_gamm_max ; ++i) {
183 ii = double(i);
184 ii2a = ii+ii+a;
185 tmp = ii*(b-ii)*x/((ii2a-1)*ii2a);
186 a0 = (a1+a0*tmp);
187 b0 = (b1+b0*tmp);
188 ii += a;
189 tmp = -ii*(ii+b)*x/((ii2a+1.)*ii2a);
190 a1 = a0+tmp*a1;
191 b1 = b0+tmp*b1;
192 mask = abs(b1-a1*rold) < abs(b1)*_gamm_err;
193 list = where(mask);
194 if (numberof(list)) r(ndx(list)) = b1(list)/a1(list);
195 list = where(!mask);
196 if (!numberof(list)) break;
197 ndx = ndx(list);
198 a = a(list);
199 b = b(list);
200 x = x(list);
201 ren = 1./a1(list);
202 a0 = a0(list)*ren;
203 b0 = b0(list)*ren;
204 b1 = rold = b1(list)*ren;
205 a1 = 1.;
206 }
207 if (numberof(list))
208 error, "global variable _gamm_err or _gamm_max too small";
209 if (numberof(r) == 1) r = r(1);
210 return r;
211 }
212