redoing dev style to be more test centric, from lessons learned with lisp-matrix.
[CommonLispStat.git] / lib / betabase.c
blob7e5e6d4a4cfbada35886e1cefd428e66a03a3b7e
1 #include "xmath.h"
3 #define TRUE 1
4 #define FALSE 0
6 /* external function */
7 extern double macheps(), gamma();
9 /* forward declarations */
10 static double logbeta(), betai(), xinbta();
12 void
13 betabase(double *x, double *a, double *b, long *gia, long *gib, double *cdf)
15 *cdf = betai(*x, *a, *b);
18 double ppbeta(double p, double a, double b, int *ifault)
20 double lbeta;
22 lbeta = gamma(a) + gamma(b) - gamma(a + b);
23 return(xinbta(&a, &b, &lbeta, &p, ifault));
27 Static routines
30 static double logbeta(double p, double q)
32 return(gamma(p) + gamma(q) - gamma(p + q));
35 #define Min(x,y) (((x) < (y)) ? (x) : (y))
36 #define Max(x,y) (((x) > (y)) ? (x) : (y))
38 static double betai(double x, double pin, double qin)
40 /* Translated from FORTRAN
41 july 1977 edition. w. fullerton, c3, los alamos scientific lab.
42 based on bosten and battiste, remark on algorithm 179, comm. acm,
43 v 17, p 153, (1974).
45 input arguments --
46 x upper limit of integration. x must be in (0,1) inclusive.
47 p first beta distribution parameter. p must be gt 0.0.
48 q second beta distribution parameter. q must be gt 0.0.
49 betai the incomplete beta function ratio is the probability that a
50 random variable from a beta distribution having parameters
51 p and q will be less than or equal to x.
53 double c, finsum, p, ps, q, term, xb, xi, y, dbetai, p1;
54 static double eps = 0.0, alneps = 0.0, sml = 0.0, alnsml = 0.0;
55 long i, n, ib;
57 /* I'm not sure these tolerances are optimal */
58 if (eps == 0.0) {
59 eps = macheps();
60 alneps = log(eps);
61 sml = macheps();
62 alnsml = log(sml);
65 y = x;
66 p = pin;
67 q = qin;
68 if (q > p || x >= 0.8)
69 if (x >= 0.2) {
70 y = 1.0 - y;
71 p = qin;
72 q = pin;
75 if ((p + q) * y / (p + 1.0) < eps) {
76 dbetai = 0.0;
77 xb = p * log(Max(y, sml)) - log(p) - logbeta(p, q);
78 if (xb > alnsml && y != 0.0) dbetai = exp(xb);
79 if (y != x || p != pin) dbetai = 1.0 - dbetai;
81 else {
83 * evaluate the infinite sum first. term will equal
84 * y**p/beta(ps,p) * (1.-ps)-sub-i * y**i / fac(i) .
86 ps = q - floor(q);
87 if (ps == 0.0) ps = 1.0;
88 xb = p * log(y) - logbeta(ps, p) - log(p);
89 dbetai = 0.0;
90 if (xb >= alnsml) {
92 dbetai = exp(xb);
93 term = dbetai * p;
94 if (ps != 1.0) {
95 n = Max(alneps / log(y), 4.0);
96 for (i = 1; i <= n; i++) {
97 xi = i;
98 term = term * (xi - ps) * y / xi;
99 dbetai = dbetai + term / (p + xi);
104 * now evaluate the finite sum, maybe.
106 if (q > 1.0) {
108 xb = p * log(y) + q * log(1.0 - y) - logbeta(p,q) - log(q);
109 ib = Max(xb / alnsml, 0.0);
110 term = exp(xb - ((float) ib) * alnsml);
111 c = 1.0 / (1.0 - y);
112 p1 = q * c / (p + q - 1.0);
114 finsum = 0.0;
115 n = q;
116 if (q == (float) n) n = n - 1;
117 for (i = 1; i <= n; i++) {
118 if (p1 <= 1.0 && term / eps <= finsum) break;
119 xi = i;
120 term = (q - xi + 1.0) * c * term / (p + q - xi);
122 if (term > 1.0) ib = ib - 1;
123 if (term > 1.0) term = term * sml;
125 if (ib==0) finsum = finsum + term;
128 dbetai = dbetai + finsum;
130 if (y != x || p != pin) dbetai = 1.0 - dbetai;
131 dbetai = Max(Min(dbetai, 1.0), 0.0);
133 return(dbetai);
137 xinbta.f -- translated by f2c and modified
139 algorithm as 109 appl. statist. (1977), vol.26, no.1
140 (replacing algorithm as 64 appl. statist. (1973), vol.22, no.3)
142 Remark AS R83 has been incorporated in this version.
144 Computes inverse of the incomplete beta function
145 ratio for given positive values of the arguments
146 p and q, alpha between zero and one.
147 log of complete beta function, beta, is assumed to be known.
149 Auxiliary function required: betai
151 SAE below is the most negative decimal exponent which does not
152 cause an underflow; a value of -308 or thereabouts will often be
155 static double xinbta(double *p, double *q, double *beta, double *alpha, int *ifault)
157 /* Initialized data */
158 static double sae = -30.0; /* this should be sufficient */
159 static double zero = 0.0;
160 static double one = 1.0;
161 static double two = 2.0;
162 static double three = 3.0;
163 static double four = 4.0;
164 static double five = 5.0;
165 static double six = 6.0;
167 /* System generated locals */
168 double ret_val, d_1, d_2;
170 /* Local variables */
171 static int indx;
172 static double prev, a, g, h, r, s, t, w, y, yprev, pp, qq;
173 static double sq, tx, adj, acu;
174 static long iex;
175 static double fpu, xin;
177 /* Define accuracy and initialise. */
178 fpu = sae * 10.;
179 ret_val = *alpha;
181 /* test for admissibility of parameters */
182 *ifault = 1;
183 if (*p <= zero || *q <= zero) return ret_val;
184 *ifault = 2;
185 if (*alpha < zero || *alpha > one) return ret_val;
186 *ifault = 0;
187 if (*alpha == zero || *alpha == one) return ret_val;
189 /* change tail if necessary */
190 if (*alpha <= .5) {
191 a = *alpha;
192 pp = *p;
193 qq = *q;
194 indx = FALSE;
196 else {
197 a = one - *alpha;
198 pp = *q;
199 qq = *p;
200 indx = TRUE;
203 /* calculate the initial approximation */
204 r = sqrt(-log(a * a));
205 y = r - (r * .27061 + 2.30753) / (one + (r * .04481 + .99229) * r);
206 if (pp > one && qq > one) {
207 r = (y * y - three) / six;
208 s = one / (pp + pp - one);
209 t = one / (qq + qq - one);
210 h = two / (s + t);
211 d_1 = y * sqrt(h + r) / h;
212 d_2 = (t - s) * (r + five / six - two / (three * h));
213 w = d_1 - d_2;
214 ret_val = pp / (pp + qq * exp(w + w));
216 else {
217 r = qq + qq;
218 t = one / (qq * 9.);
219 /* Computing 3rd power */
220 d_1 = one - t + y * sqrt(t);
221 t = r * (d_1 * d_1 * d_1);
222 if (t <= zero) {
223 ret_val = one - exp((log((one - a) * qq) + *beta) / qq);
225 else {
226 t = (four * pp + r - two) / t;
227 if (t <= one) ret_val = exp((log(a * pp) + *beta) / pp);
228 else ret_val = one - two / (t + one);
234 solve for x by a modified newton-raphson method, using the function betai
236 r = one - pp;
237 t = one - qq;
238 yprev = zero;
239 sq = one;
240 prev = one;
241 if (ret_val < 1e-4) ret_val = 1e-4;
242 if (ret_val > .9999) ret_val = .9999;
243 /* Computing MAX, two 2nd powers */
244 d_1 = -5.0 / (pp * pp) - 1.0 / (a * a) - 13.0;
245 iex = (sae > d_1) ? sae : d_1;
246 acu = pow(10.0, (double) iex);
247 do {
248 y = betai(ret_val, pp, qq);
249 if (*ifault != 0) {
250 *ifault = 3;
251 return ret_val;
253 xin = ret_val;
254 y = (y - a) * exp(*beta + r * log(xin) + t * log(one - xin));
255 if (y * yprev <= zero) {
256 prev = (sq > fpu) ? sq : fpu;
258 g = one;
259 do {
260 adj = g * y;
261 sq = adj * adj;
262 if (sq < prev) {
263 tx = ret_val - adj;
264 if (tx >= zero && tx <= one) {
265 if (prev <= acu || y * y <= acu) {
266 if (indx) ret_val = one - ret_val;
267 return ret_val;
269 if (tx != zero && tx != one) break;
272 g /= three;
273 } while (TRUE);
274 if (tx == ret_val) {
275 if (indx) ret_val = one - ret_val;
276 return ret_val;
278 ret_val = tx;
279 yprev = y;
280 } while (TRUE);
281 return ret_val;