6 /* external function */
7 extern double macheps(), gamma();
9 /* forward declarations */
10 static double logbeta(), betai(), xinbta();
13 betabase(double *x
, double *a
, double *b
, int *gia
, int *gib
, double *cdf
)
15 *cdf
= betai(*x
, *a
, *b
);
18 double ppbeta(double p
, double a
, double b
, int *ifault
)
22 lbeta
= gamma(a
) + gamma(b
) - gamma(a
+ b
);
23 return(xinbta(&a
, &b
, &lbeta
, &p
, ifault
));
30 static double logbeta(p
, q
)
33 return(gamma(p
) + gamma(q
) - gamma(p
+ q
));
36 #define Min(x,y) (((x) < (y)) ? (x) : (y))
37 #define Max(x,y) (((x) > (y)) ? (x) : (y))
39 static double betai(x
, pin
, qin
)
42 /* Translated from FORTRAN
43 july 1977 edition. w. fullerton, c3, los alamos scientific lab.
44 based on bosten and battiste, remark on algorithm 179, comm. acm,
48 x upper limit of integration. x must be in (0,1) inclusive.
49 p first beta distribution parameter. p must be gt 0.0.
50 q second beta distribution parameter. q must be gt 0.0.
51 betai the incomplete beta function ratio is the probability that a
52 random variable from a beta distribution having parameters
53 p and q will be less than or equal to x.
55 double c
, finsum
, p
, ps
, q
, term
, xb
, xi
, y
, dbetai
, p1
;
56 static double eps
= 0.0, alneps
= 0.0, sml
= 0.0, alnsml
= 0.0;
59 /* I'm not sure these tolerances are optimal */
70 if (q
> p
|| x
>= 0.8)
77 if ((p
+ q
) * y
/ (p
+ 1.0) < eps
) {
79 xb
= p
* log(Max(y
, sml
)) - log(p
) - logbeta(p
, q
);
80 if (xb
> alnsml
&& y
!= 0.0) dbetai
= exp(xb
);
81 if (y
!= x
|| p
!= pin
) dbetai
= 1.0 - dbetai
;
85 * evaluate the infinite sum first. term will equal
86 * y**p/beta(ps,p) * (1.-ps)-sub-i * y**i / fac(i) .
89 if (ps
== 0.0) ps
= 1.0;
90 xb
= p
* log(y
) - logbeta(ps
, p
) - log(p
);
97 n
= Max(alneps
/ log(y
), 4.0);
98 for (i
= 1; i
<= n
; i
++) {
100 term
= term
* (xi
- ps
) * y
/ xi
;
101 dbetai
= dbetai
+ term
/ (p
+ xi
);
106 * now evaluate the finite sum, maybe.
110 xb
= p
* log(y
) + q
* log(1.0 - y
) - logbeta(p
,q
) - log(q
);
111 ib
= Max(xb
/ alnsml
, 0.0);
112 term
= exp(xb
- ((float) ib
) * alnsml
);
114 p1
= q
* c
/ (p
+ q
- 1.0);
118 if (q
== (float) n
) n
= n
- 1;
119 for (i
= 1; i
<= n
; i
++) {
120 if (p1
<= 1.0 && term
/ eps
<= finsum
) break;
122 term
= (q
- xi
+ 1.0) * c
* term
/ (p
+ q
- xi
);
124 if (term
> 1.0) ib
= ib
- 1;
125 if (term
> 1.0) term
= term
* sml
;
127 if (ib
==0) finsum
= finsum
+ term
;
130 dbetai
= dbetai
+ finsum
;
132 if (y
!= x
|| p
!= pin
) dbetai
= 1.0 - dbetai
;
133 dbetai
= Max(Min(dbetai
, 1.0), 0.0);
139 xinbta.f -- translated by f2c and modified
141 algorithm as 109 appl. statist. (1977), vol.26, no.1
142 (replacing algorithm as 64 appl. statist. (1973), vol.22, no.3)
144 Remark AS R83 has been incorporated in this version.
146 Computes inverse of the incomplete beta function
147 ratio for given positive values of the arguments
148 p and q, alpha between zero and one.
149 log of complete beta function, beta, is assumed to be known.
151 Auxiliary function required: betai
153 SAE below is the most negative decimal exponent which does not
154 cause an underflow; a value of -308 or thereabouts will often be
157 static double xinbta(p
, q
, beta
, alpha
, ifault
)
158 double *p
, *q
, *beta
, *alpha
;
161 /* Initialized data */
162 static double sae
= -30.0; /* this should be sufficient */
163 static double zero
= 0.0;
164 static double one
= 1.0;
165 static double two
= 2.0;
166 static double three
= 3.0;
167 static double four
= 4.0;
168 static double five
= 5.0;
169 static double six
= 6.0;
171 /* System generated locals */
172 double ret_val
, d_1
, d_2
;
174 /* Local variables */
176 static double prev
, a
, g
, h
, r
, s
, t
, w
, y
, yprev
, pp
, qq
;
177 static double sq
, tx
, adj
, acu
;
179 static double fpu
, xin
;
181 /* Define accuracy and initialise. */
185 /* test for admissibility of parameters */
187 if (*p
<= zero
|| *q
<= zero
) return ret_val
;
189 if (*alpha
< zero
|| *alpha
> one
) return ret_val
;
191 if (*alpha
== zero
|| *alpha
== one
) return ret_val
;
193 /* change tail if necessary */
207 /* calculate the initial approximation */
208 r
= sqrt(-log(a
* a
));
209 y
= r
- (r
* .27061 + 2.30753) / (one
+ (r
* .04481 + .99229) * r
);
210 if (pp
> one
&& qq
> one
) {
211 r
= (y
* y
- three
) / six
;
212 s
= one
/ (pp
+ pp
- one
);
213 t
= one
/ (qq
+ qq
- one
);
215 d_1
= y
* sqrt(h
+ r
) / h
;
216 d_2
= (t
- s
) * (r
+ five
/ six
- two
/ (three
* h
));
218 ret_val
= pp
/ (pp
+ qq
* exp(w
+ w
));
223 /* Computing 3rd power */
224 d_1
= one
- t
+ y
* sqrt(t
);
225 t
= r
* (d_1
* d_1
* d_1
);
227 ret_val
= one
- exp((log((one
- a
) * qq
) + *beta
) / qq
);
230 t
= (four
* pp
+ r
- two
) / t
;
231 if (t
<= one
) ret_val
= exp((log(a
* pp
) + *beta
) / pp
);
232 else ret_val
= one
- two
/ (t
+ one
);
238 solve for x by a modified newton-raphson method, using the function betai
245 if (ret_val
< 1e-4) ret_val
= 1e-4;
246 if (ret_val
> .9999) ret_val
= .9999;
247 /* Computing MAX, two 2nd powers */
248 d_1
= -5.0 / (pp
* pp
) - 1.0 / (a
* a
) - 13.0;
249 iex
= (sae
> d_1
) ? sae
: d_1
;
250 acu
= pow(10.0, (double) iex
);
252 y
= betai(ret_val
, pp
, qq
);
258 y
= (y
- a
) * exp(*beta
+ r
* log(xin
) + t
* log(one
- xin
));
259 if (y
* yprev
<= zero
) {
260 prev
= (sq
> fpu
) ? sq
: fpu
;
268 if (tx
>= zero
&& tx
<= one
) {
269 if (prev
<= acu
|| y
* y
<= acu
) {
270 if (indx
) ret_val
= one
- ret_val
;
273 if (tx
!= zero
&& tx
!= one
) break;
279 if (indx
) ret_val
= one
- ret_val
;