3 * Copyright (C) 2004 Jason H. Stover.
5 * This program is free software; you can redistribute it and/or modify
6 * it under the terms of the GNU General Public License as published by
7 * the Free Software Foundation; either version 2 of the License, or (at
8 * your option) any later version.
10 * This program is distributed in the hope that it will be useful, but
11 * WITHOUT ANY WARRANTY; without even the implied warranty of
12 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
13 * General Public License for more details.
15 * You should have received a copy of the GNU General Public License
16 * along with this program; if not, write to the Free Software
17 * Foundation, Inc., 59 Temple Place, Suite 330, Boston, MA 02111-1307, USA.
21 * Invert the Beta distribution.
25 * Roger W. Abernathy and Robert P. Smith. "Applying Series Expansion
26 * to the Inverse Beta Distribution to Find Percentiles of the F-Distribution,"
27 * ACM Transactions on Mathematical Software, volume 19, number 4, December 1993,
30 * G.W. Hill and A.W. Davis. "Generalized asymptotic expansions of a
31 * Cornish-Fisher type," Annals of Mathematical Statistics, volume 39, number 8,
32 * August 1968, pages 1264-1273.
36 #include <gsl/gsl_math.h>
37 #include <gsl/gsl_errno.h>
38 #include <gsl/gsl_sf_gamma.h>
39 #include <gsl/gsl_cdf.h>
40 #include <gsl/gsl_randist.h>
41 #include "gsl-extras.h"
43 #define BETAINV_INIT_ERR .01
44 #define BETADISTINV_N_TERMS 3
45 #define BETADISTINV_MAXITER 20
48 s_bisect (double x, double y)
50 double result = GSL_MIN(x,y) + fabs(x - y) / 2.0;
54 new_guess_P ( double old_guess, double x, double y,
55 double prob, double a, double b)
61 p_hat = gsl_cdf_beta_P(old_guess, a, b);
64 end_point = GSL_MAX(x,y);
66 else if ( p_hat > prob )
68 end_point = GSL_MIN(x,y);
72 end_point = old_guess;
74 result = s_bisect(old_guess, end_point);
80 new_guess_Q ( double old_guess, double x, double y,
81 double prob, double a, double b)
87 q_hat = gsl_cdf_beta_Q(old_guess, a, b);
90 end_point = GSL_MAX(x,y);
92 else if ( q_hat < prob )
94 end_point = GSL_MIN(x,y);
98 end_point = old_guess;
100 result = s_bisect(old_guess, end_point);
106 * The get_corn_fish_* functions below return the first
107 * three terms of the Cornish-Fisher expansion
108 * without recursion. The recursive functions
109 * make the code more legible when higher order coefficients
110 * are used, but terms beyond the cubic do not
114 * Linear coefficient for the
115 * Cornish-Fisher expansion.
118 get_corn_fish_lin (const double x, const double a, const double b)
122 result = gsl_ran_beta_pdf (x, a, b);
125 result = 1.0 / result;
129 result = GSL_DBL_MAX;
135 * Quadratic coefficient for the
136 * Cornish-Fisher expansion.
139 get_corn_fish_quad (const double x, const double a, const double b)
148 gam_ab = gsl_sf_lngamma(a + b);
149 gam_a = gsl_sf_lngamma (a);
150 gam_b = gsl_sf_lngamma (b);
151 num = exp(2 * (gam_a + gam_b - gam_ab)) * (1 - a + x * (b + a - 2));
152 den = 2.0 * pow ( x, 2*a - 1 ) * pow ( 1 - x, 2 * b - 1 );
153 if ( fabs(den) > 0.0)
165 * The cubic term for the Cornish-Fisher expansion.
166 * Theoretically, use of this term should give a better approximation,
167 * but in practice inclusion of the cubic term worsens the
168 * iterative procedure in gsl_cdf_beta_Pinv and gsl_cdf_beta_Qinv
169 * for extreme values of p, a or b.
173 get_corn_fish_cube (const double x, const double a, const double b)
176 double am1 = a - 1.0;
177 double am2 = a - 2.0;
178 double apbm2 = a+b-2.0;
179 double apbm3 = a+b-3.0;
180 double apbm4 = a+b-4.0;
181 double ab1ab2 = am1 * am2;
186 num = (am1 - x * apbm2) * (am1 - x * apbm2);
187 tmp = ab1ab2 - x * (apbm2 * ( ab1ab2 * apbm4 + 1) + x * apbm2 * apbm3);
189 tmp = gsl_ran_beta_pdf(x,a,b);
190 den = 2.0 * x * x * (1 - x) * (1 - x) * pow(tmp,3.0);
191 if ( fabs(den) > 0.0)
204 * The Cornish-Fisher coefficients can be defined recursivley,
205 * starting with the nth derivative of s_psi = -f'(x)/f(x),
206 * where f is the beta density.
208 * The section below was commented out since
209 * the recursive generation of the coeficients did
210 * not improve the accuracy of the directly coded
211 * the first three coefficients.
215 s_d_psi (double x, double a, double b, int n)
218 double np1 = (double) n + 1;
221 double bm1 = b - 1.0;
222 double am1 = a - 1.0;
225 asgn = (n % 2) ? 1.0:-1.0;
226 bsgn = (n % 2) ? -1.0:1.0;
227 result = gsl_sf_gamma(np1) * ((bsgn * bm1 / (pow(mx, np1)))
228 + (asgn * am1 / (pow(x,np1))));
232 * nth derivative of a coefficient with respect
236 get_d_coeff ( double x, double a,
237 double b, double n, double k)
248 result = s_d_psi(x, a, b, k);
253 for (i = 0.0; i < (k+1); i++)
255 k_fac = gsl_sf_lngamma(k+1.0);
256 i_fac = gsl_sf_lngamma(i+1.0);
257 kmi_fac = gsl_sf_lngamma(k-i+1.0);
259 result += exp(k_fac - i_fac - kmi_fac)
260 * get_d_coeff( x, a, b, 2.0, i)
261 * get_d_coeff( x, a, b, (n - 1.0), (k - i));
263 result += get_d_coeff ( x, a, b, (n-1.0), (k+1.0));
269 * Cornish-Fisher coefficient.
272 get_corn_fish (double c, double x,
273 double a, double b, double n)
285 result = s_d_psi(x, a, b, 0);
289 dc = get_d_coeff(x, a, b, (n-1.0), 1.0);
290 c_prev = get_corn_fish(c, x, a, b, (n-1.0));
291 result = (n-1) * s_d_psi(x,a,b,0) * c_prev + dc;
298 gslextras_cdf_beta_Pinv ( const double p, const double a, const double b)
326 GSLEXTRAS_CDF_ERROR("p < 0", GSL_EDOM);
330 GSLEXTRAS_CDF_ERROR("p > 1",GSL_EDOM);
334 GSLEXTRAS_CDF_ERROR ("a < 0", GSL_EDOM );
338 GSLEXTRAS_CDF_ERROR ( "b < 0", GSL_EDOM );
349 if (p > (1.0 - GSL_DBL_EPSILON))
352 * When p is close to 1.0, the bisection
353 * works better with gsl_cdf_Q.
355 state = gslextras_cdf_beta_Qinv ( p, a, b);
356 result = 1.0 - state;
359 if (p < GSL_DBL_EPSILON )
362 * Start at a small value and rise until
363 * we are above the correct result. This
364 * avoids overflow. When p is very close to
365 * 0, an initial state value of a/(a+b) will
366 * cause the interpolating polynomial
370 beta_result = gsl_cdf_beta_P (upper, a, b);
371 while (beta_result < p)
375 beta_result = gsl_cdf_beta_P (upper, a, b);
377 state = (lower + upper) / 2.0;
382 * First guess is the expected value.
387 beta_result = gsl_cdf_beta_P (state, a, b);
389 err = beta_result - p;
392 while ( relerr > BETAINV_INIT_ERR && n_iter < 100)
394 tmp = new_guess_P ( state, lower, upper,
396 lower = ( tmp < state ) ? lower:state;
397 upper = ( tmp < state ) ? state:upper;
399 beta_result = gsl_cdf_beta_P (state, a, b);
400 err = p - beta_result;
408 * Use a second order Lagrange interpolating
409 * polynomial to get closer before switching to
410 * the iterative method.
412 p0 = gsl_cdf_beta_P (lower, a, b);
413 p1 = gsl_cdf_beta_P (state, a, b);
414 p2 = gsl_cdf_beta_P (upper, a, b);
415 if( p0 < p1 && p1 < p2)
417 frac1 = (p - p2) / (p0 - p1);
418 frac2 = (p - p1) / (p0 - p2);
419 frac3 = (p - p0) / (p1 - p2);
420 frac4 = (p - p0) * (p - p1) / ((p2 - p0) * (p2 - p1));
421 state = frac1 * (frac2 * lower - frac3 * state)
424 beta_result = gsl_cdf_beta_P( state, a, b);
425 err = p - beta_result;
428 if (relerr < min_err)
438 * Newton-type iteration using the terms from the
439 * Cornish-Fisher expansion. If only the first term
440 * of the expansion is used, this is Newton's method.
442 while ( relerr > GSL_DBL_EPSILON && n_iter < BETADISTINV_MAXITER)
445 c1 = get_corn_fish_lin (state, a, b);
446 c2 = get_corn_fish_quad (state, a, b);
448 * The cubic term does not help, and can can
449 * harm the approximation for extreme values of
453 c3 = get_corn_fish_cube (state, a, b);
454 state += err * (c1 + (err / 2.0 ) * (c2 + c3 * err / 3.0));
456 state += err * (c1 + (c2 * err / 2.0 ));
458 * The section below which is commented out uses
459 * a recursive function to get the coefficients.
460 * The recursion makes coding higher-order terms
461 * easier, but did not improve the result beyond
462 * the use of three terms. Since explicitly coding
463 * those three terms in the get_corn_fish_* functions
464 * was not difficult, the recursion was abandoned.
468 for(i = 1.0; i < BETADISTINV_N_TERMS; i += 1.0)
471 coeff = get_corn_fish (coeff, prior_state, a, b, i);
472 state += coeff * pow(err, i) /
473 (i_fac * pow (gsl_ran_beta_pdf(prior_state,a,b), i));
476 beta_result = gsl_cdf_beta_P ( state, a, b );
477 err = p - beta_result;
480 if (relerr < min_err)
491 gslextras_cdf_beta_Qinv (double q, double a, double b)
519 GSLEXTRAS_CDF_ERROR("q < 0", GSL_EDOM);
523 GSLEXTRAS_CDF_ERROR("q > 1",GSL_EDOM);
527 GSLEXTRAS_CDF_ERROR ("a < 0", GSL_EDOM );
531 GSLEXTRAS_CDF_ERROR ( "b < 0", GSL_EDOM );
542 if ( q < GSL_DBL_EPSILON )
545 * When q is close to 0, the bisection
546 * and interpolation done in the rest of
547 * this routine will not give the correct
548 * value within double precision, so
549 * gsl_cdf_beta_Qinv is called instead.
551 state = gslextras_cdf_beta_Pinv ( q, a, b);
552 result = 1.0 - state;
555 if ( q > 1.0 - GSL_DBL_EPSILON )
558 * Make the initial guess close to 0.0.
561 beta_result = gsl_cdf_beta_Q ( upper, a, b);
562 while (beta_result > q )
566 beta_result = gsl_cdf_beta_Q ( upper, a, b);
568 state = (upper + lower) / 2.0;
572 /* Bisection to get an initial approximation.
573 * First guess is the expected value.
579 beta_result = gsl_cdf_beta_Q (state, a, b);
580 err = beta_result - q;
583 while ( relerr > BETAINV_INIT_ERR && n_iter < 100)
586 tmp = new_guess_Q ( state, lower, upper,
588 lower = ( tmp < state ) ? lower:state;
589 upper = ( tmp < state ) ? state:upper;
591 beta_result = gsl_cdf_beta_Q (state, a, b);
592 err = q - beta_result;
600 * Use a second order Lagrange interpolating
601 * polynomial to get closer before switching to
602 * the iterative method.
604 p0 = gsl_cdf_beta_Q (lower, a, b);
605 p1 = gsl_cdf_beta_Q (state, a, b);
606 p2 = gsl_cdf_beta_Q (upper, a, b);
607 if(p0 > p1 && p1 > p2)
609 frac1 = (q - p2) / (p0 - p1);
610 frac2 = (q - p1) / (p0 - p2);
611 frac3 = (q - p0) / (p1 - p2);
612 frac4 = (q - p0) * (q - p1) / ((p2 - p0) * (p2 - p1));
613 state = frac1 * (frac2 * lower - frac3 * state)
615 beta_result = gsl_cdf_beta_Q( state, a, b);
616 err = beta_result - q;
619 if (relerr < min_err)
627 * Iteration using the terms from the
628 * Cornish-Fisher expansion. If only the first term
629 * of the expansion is used, this is Newton's method.
633 while ( relerr > GSL_DBL_EPSILON && n_iter < BETADISTINV_MAXITER)
636 c1 = get_corn_fish_lin (state, a, b);
637 c2 = get_corn_fish_quad (state, a, b);
639 * The cubic term does not help, and can harm
640 * the approximation for extreme values of p, a and b.
643 c3 = get_corn_fish_cube (state, a, b);
644 state += err * (c1 + (err / 2.0 ) * (c2 + c3 * err / 3.0));
646 state += err * (c1 + (c2 * err / 2.0 ));
647 beta_result = gsl_cdf_beta_Q ( state, a, b );
648 err = beta_result - q;
651 if (relerr < min_err)