3 * Copyright (C) 2004 Free Software Foundation, Inc.
4 * Written by Jason H. Stover.
6 * This program is free software; you can redistribute it and/or modify
7 * it under the terms of the GNU General Public License as published by
8 * the Free Software Foundation; either version 2 of the License, or (at
9 * your option) any later version.
11 * This program is distributed in the hope that it will be useful, but
12 * WITHOUT ANY WARRANTY; without even the implied warranty of
13 * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the GNU
14 * General Public License for more details.
16 * You should have received a copy of the GNU General Public License
17 * along with this program; if not, write to the Free Software
18 * Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
22 * Invert the Beta distribution.
26 * Roger W. Abernathy and Robert P. Smith. "Applying Series Expansion
27 * to the Inverse Beta Distribution to Find Percentiles of the F-Distribution,"
28 * ACM Transactions on Mathematical Software, volume 19, number 4, December 1993,
31 * G.W. Hill and A.W. Davis. "Generalized asymptotic expansions of a
32 * Cornish-Fisher type," Annals of Mathematical Statistics, volume 39, number 8,
33 * August 1968, pages 1264-1273.
38 #include <gsl/gsl_math.h>
39 #include <gsl/gsl_errno.h>
40 #include <gsl/gsl_sf_gamma.h>
41 #include <gsl/gsl_cdf.h>
42 #include <gsl/gsl_randist.h>
43 #include "gsl-extras.h"
45 #define BETAINV_INIT_ERR .001
46 #define BETADISTINV_N_TERMS 3
47 #define BETADISTINV_MAXITER 20
50 s_bisect (double x, double y)
52 double result = GSL_MIN(x,y) + fabs(x - y) / 2.0;
56 new_guess_P ( double old_guess, double x, double y,
57 double prob, double a, double b)
63 p_hat = gsl_cdf_beta_P(old_guess, a, b);
66 end_point = GSL_MAX(x,y);
68 else if ( p_hat > prob )
70 end_point = GSL_MIN(x,y);
74 end_point = old_guess;
76 result = s_bisect(old_guess, end_point);
82 new_guess_Q ( double old_guess, double x, double y,
83 double prob, double a, double b)
89 q_hat = gsl_cdf_beta_Q(old_guess, a, b);
92 end_point = GSL_MAX(x,y);
94 else if ( q_hat < prob )
96 end_point = GSL_MIN(x,y);
100 end_point = old_guess;
102 result = s_bisect(old_guess, end_point);
108 * The get_corn_fish_* functions below return the first
109 * three terms of the Cornish-Fisher expansion
110 * without recursion. The recursive functions
111 * make the code more legible when higher order coefficients
112 * are used, but terms beyond the cubic do not
116 * Linear coefficient for the
117 * Cornish-Fisher expansion.
120 get_corn_fish_lin (const double x, const double a, const double b)
124 result = gsl_ran_beta_pdf (x, a, b);
127 result = 1.0 / result;
131 result = GSL_DBL_MAX;
137 * Quadratic coefficient for the
138 * Cornish-Fisher expansion.
141 get_corn_fish_quad (const double x, const double a, const double b)
150 gam_ab = gsl_sf_lngamma(a + b);
151 gam_a = gsl_sf_lngamma (a);
152 gam_b = gsl_sf_lngamma (b);
153 num = exp(2 * (gam_a + gam_b - gam_ab)) * (1 - a + x * (b + a - 2));
154 den = 2.0 * pow ( x, 2*a - 1 ) * pow ( 1 - x, 2 * b - 1 );
155 if ( fabs(den) > 0.0)
167 * The cubic term for the Cornish-Fisher expansion.
168 * Theoretically, use of this term should give a better approximation,
169 * but in practice inclusion of the cubic term worsens the
170 * iterative procedure in gsl_cdf_beta_Pinv and gsl_cdf_beta_Qinv
171 * for extreme values of p, a or b.
175 get_corn_fish_cube (const double x, const double a, const double b)
178 double am1 = a - 1.0;
179 double am2 = a - 2.0;
180 double apbm2 = a+b-2.0;
181 double apbm3 = a+b-3.0;
182 double apbm4 = a+b-4.0;
183 double ab1ab2 = am1 * am2;
188 num = (am1 - x * apbm2) * (am1 - x * apbm2);
189 tmp = ab1ab2 - x * (apbm2 * ( ab1ab2 * apbm4 + 1) + x * apbm2 * apbm3);
191 tmp = gsl_ran_beta_pdf(x,a,b);
192 den = 2.0 * x * x * (1 - x) * (1 - x) * pow(tmp,3.0);
193 if ( fabs(den) > 0.0)
206 * The Cornish-Fisher coefficients can be defined recursivley,
207 * starting with the nth derivative of s_psi = -f'(x)/f(x),
208 * where f is the beta density.
210 * The section below was commented out since
211 * the recursive generation of the coeficients did
212 * not improve the accuracy of the directly coded
213 * the first three coefficients.
217 s_d_psi (double x, double a, double b, int n)
220 double np1 = (double) n + 1;
223 double bm1 = b - 1.0;
224 double am1 = a - 1.0;
227 asgn = (n % 2) ? 1.0:-1.0;
228 bsgn = (n % 2) ? -1.0:1.0;
229 result = gsl_sf_gamma(np1) * ((bsgn * bm1 / (pow(mx, np1)))
230 + (asgn * am1 / (pow(x,np1))));
234 * nth derivative of a coefficient with respect
238 get_d_coeff ( double x, double a,
239 double b, double n, double k)
250 result = s_d_psi(x, a, b, k);
255 for (i = 0.0; i < (k+1); i++)
257 k_fac = gsl_sf_lngamma(k+1.0);
258 i_fac = gsl_sf_lngamma(i+1.0);
259 kmi_fac = gsl_sf_lngamma(k-i+1.0);
261 result += exp(k_fac - i_fac - kmi_fac)
262 * get_d_coeff( x, a, b, 2.0, i)
263 * get_d_coeff( x, a, b, (n - 1.0), (k - i));
265 result += get_d_coeff ( x, a, b, (n-1.0), (k+1.0));
271 * Cornish-Fisher coefficient.
274 get_corn_fish (double c, double x,
275 double a, double b, double n)
287 result = s_d_psi(x, a, b, 0);
291 dc = get_d_coeff(x, a, b, (n-1.0), 1.0);
292 c_prev = get_corn_fish(c, x, a, b, (n-1.0));
293 result = (n-1) * s_d_psi(x,a,b,0) * c_prev + dc;
300 gslextras_cdf_beta_Pinv ( const double p, const double a, const double b)
328 GSLEXTRAS_CDF_ERROR("p < 0", GSL_EDOM);
332 GSLEXTRAS_CDF_ERROR("p > 1",GSL_EDOM);
336 GSLEXTRAS_CDF_ERROR ("a < 0", GSL_EDOM );
340 GSLEXTRAS_CDF_ERROR ( "b < 0", GSL_EDOM );
351 if (p > (1.0 - GSL_DBL_EPSILON))
354 * When p is close to 1.0, the bisection
355 * works better with gsl_cdf_Q.
357 state = gslextras_cdf_beta_Qinv ( p, a, b);
358 result = 1.0 - state;
361 if (p < GSL_DBL_EPSILON )
364 * Start at a small value and rise until
365 * we are above the correct result. This
366 * avoids overflow. When p is very close to
367 * 0, an initial state value of a/(a+b) will
368 * cause the interpolating polynomial
372 beta_result = gsl_cdf_beta_P (upper, a, b);
373 while (beta_result < p)
377 beta_result = gsl_cdf_beta_P (upper, a, b);
379 state = (lower + upper) / 2.0;
384 * First guess is the expected value.
389 beta_result = gsl_cdf_beta_P (state, a, b);
391 err = beta_result - p;
394 while ( relerr > BETAINV_INIT_ERR)
396 tmp = new_guess_P ( state, lower, upper,
398 lower = ( tmp < state ) ? lower:state;
399 upper = ( tmp < state ) ? state:upper;
401 beta_result = gsl_cdf_beta_P (state, a, b);
402 err = p - beta_result;
410 * Use a second order Lagrange interpolating
411 * polynomial to get closer before switching to
412 * the iterative method.
414 p0 = gsl_cdf_beta_P (lower, a, b);
415 p1 = gsl_cdf_beta_P (state, a, b);
416 p2 = gsl_cdf_beta_P (upper, a, b);
417 if( p0 < p1 && p1 < p2)
419 frac1 = (p - p2) / (p0 - p1);
420 frac2 = (p - p1) / (p0 - p2);
421 frac3 = (p - p0) / (p1 - p2);
422 frac4 = (p - p0) * (p - p1) / ((p2 - p0) * (p2 - p1));
423 state = frac1 * (frac2 * lower - frac3 * state)
426 beta_result = gsl_cdf_beta_P( state, a, b);
427 err = p - beta_result;
430 if (relerr < min_err)
438 * Lagrange polynomial failed to reduce the error.
439 * This will happen with a very skewed beta density.
440 * Undo previous steps.
443 beta_result = gsl_cdf_beta_P(state,a,b);
444 err = p - beta_result;
453 * Newton-type iteration using the terms from the
454 * Cornish-Fisher expansion. If only the first term
455 * of the expansion is used, this is Newton's method.
457 while ( relerr > GSL_DBL_EPSILON && n_iter < BETADISTINV_MAXITER)
460 c1 = get_corn_fish_lin (state, a, b);
461 c2 = get_corn_fish_quad (state, a, b);
463 * The cubic term does not help, and can can
464 * harm the approximation for extreme values of
468 c3 = get_corn_fish_cube (state, a, b);
469 state += err * (c1 + (err / 2.0 ) * (c2 + c3 * err / 3.0));
471 state += err * (c1 + (c2 * err / 2.0 ));
473 * The section below which is commented out uses
474 * a recursive function to get the coefficients.
475 * The recursion makes coding higher-order terms
476 * easier, but did not improve the result beyond
477 * the use of three terms. Since explicitly coding
478 * those three terms in the get_corn_fish_* functions
479 * was not difficult, the recursion was abandoned.
483 for(i = 1.0; i < BETADISTINV_N_TERMS; i += 1.0)
486 coeff = get_corn_fish (coeff, prior_state, a, b, i);
487 state += coeff * pow(err, i) /
488 (i_fac * pow (gsl_ran_beta_pdf(prior_state,a,b), i));
491 beta_result = gsl_cdf_beta_P ( state, a, b );
492 err = p - beta_result;
495 if (relerr < min_err)
506 gslextras_cdf_beta_Qinv (double q, double a, double b)
534 GSLEXTRAS_CDF_ERROR("q < 0", GSL_EDOM);
538 GSLEXTRAS_CDF_ERROR("q > 1",GSL_EDOM);
542 GSLEXTRAS_CDF_ERROR ("a < 0", GSL_EDOM );
546 GSLEXTRAS_CDF_ERROR ( "b < 0", GSL_EDOM );
557 if ( q < GSL_DBL_EPSILON )
560 * When q is close to 0, the bisection
561 * and interpolation done in the rest of
562 * this routine will not give the correct
563 * value within double precision, so
564 * gsl_cdf_beta_Qinv is called instead.
566 state = gslextras_cdf_beta_Pinv ( q, a, b);
567 result = 1.0 - state;
570 if ( q > 1.0 - GSL_DBL_EPSILON )
573 * Make the initial guess close to 0.0.
576 beta_result = gsl_cdf_beta_Q ( upper, a, b);
577 while (beta_result > q )
581 beta_result = gsl_cdf_beta_Q ( upper, a, b);
583 state = (upper + lower) / 2.0;
587 /* Bisection to get an initial approximation.
588 * First guess is the expected value.
594 beta_result = gsl_cdf_beta_Q (state, a, b);
595 err = beta_result - q;
598 while ( relerr > BETAINV_INIT_ERR)
601 tmp = new_guess_Q ( state, lower, upper,
603 lower = ( tmp < state ) ? lower:state;
604 upper = ( tmp < state ) ? state:upper;
606 beta_result = gsl_cdf_beta_Q (state, a, b);
607 err = q - beta_result;
615 * Use a second order Lagrange interpolating
616 * polynomial to get closer before switching to
617 * the iterative method.
619 p0 = gsl_cdf_beta_Q (lower, a, b);
620 p1 = gsl_cdf_beta_Q (state, a, b);
621 p2 = gsl_cdf_beta_Q (upper, a, b);
622 if(p0 > p1 && p1 > p2)
624 frac1 = (q - p2) / (p0 - p1);
625 frac2 = (q - p1) / (p0 - p2);
626 frac3 = (q - p0) / (p1 - p2);
627 frac4 = (q - p0) * (q - p1) / ((p2 - p0) * (p2 - p1));
628 state = frac1 * (frac2 * lower - frac3 * state)
630 beta_result = gsl_cdf_beta_Q( state, a, b);
631 err = beta_result - q;
634 if (relerr < min_err)
642 * Lagrange polynomial failed to reduce the error.
643 * This will happen with a very skewed beta density.
644 * Undo previous steps.
647 beta_result = gsl_cdf_beta_P(state,a,b);
648 err = q - beta_result;
655 * Iteration using the terms from the
656 * Cornish-Fisher expansion. If only the first term
657 * of the expansion is used, this is Newton's method.
661 while ( relerr > GSL_DBL_EPSILON && n_iter < BETADISTINV_MAXITER)
664 c1 = get_corn_fish_lin (state, a, b);
665 c2 = get_corn_fish_quad (state, a, b);
667 * The cubic term does not help, and can harm
668 * the approximation for extreme values of p, a and b.
671 c3 = get_corn_fish_cube (state, a, b);
672 state += err * (c1 + (err / 2.0 ) * (c2 + c3 * err / 3.0));
674 state += err * (c1 + (c2 * err / 2.0 ));
675 beta_result = gsl_cdf_beta_Q ( state, a, b );
676 err = beta_result - q;
679 if (relerr < min_err)