At Fri, 18 Apr 2008 10:24:57 -0700, Bret Beck wrote: > > While testing my use of gsl_cdf_chisq_Pinv I ran into the following > error message. > > gsl: gammainv.c:105: ERROR: inverse failed to converge > Default GSL error handler invoked. > Aborted
This is now fixed in the repository at https://savannah.gnu.org/projects/gsl Thanks for the bug report. -- Brian Gough GNU Scientific Library - http://www.gnu.org/software/gsl/ commit 6d5d363d069aa108b351dff86872f7e10137f496 Author: Brian Gough <[EMAIL PROTECTED]> Date: Tue Apr 29 14:39:15 2008 +0100 fix for bug #23101 diff --git a/NEWS b/NEWS index 9f46542..e8d918e 100644 --- a/NEWS +++ b/NEWS @@ -1,3 +1,7 @@ +** Fixed a problem with convergence for inverse gamma and chisq + distribitions, gsl_cdf_gamma_{P,Q}inv and gsl_cdf_chisq_{P,Q}inv + (bug #23101). + ** Improved the handling of constant regions in Vegas by eliminating spurious excess precision when computing box variances. diff --git a/cdf/ChangeLog b/cdf/ChangeLog index a1ccdac..4caa6b6 100644 --- a/cdf/ChangeLog +++ b/cdf/ChangeLog @@ -1,3 +1,8 @@ +2008-04-29 Brian Gough <[EMAIL PROTECTED]> + + * gammainv.c (gsl_cdf_gamma_Pinv, gsl_cdf_gamma_Qinv): restrict + the range of the gaussian approximation + 2008-02-20 Brian Gough <[EMAIL PROTECTED]> * beta_inc.c (beta_inc_AXPY): add some handling for large diff --git a/cdf/gammainv.c b/cdf/gammainv.c index b88117b..767a4c0 100644 --- a/cdf/gammainv.c +++ b/cdf/gammainv.c @@ -42,7 +42,13 @@ gsl_cdf_gamma_Pinv (const double P, const double a, const double b) /* Consider, small, large and intermediate cases separately. The boundaries at 0.05 and 0.95 have not been optimised, but seem ok - for an initial approximation. */ + for an initial approximation. + + BJG: These approximations aren't really valid, the relevant + criterion is P*gamma(a+1) < 1. Need to rework these routines and + use a single bisection style solver for all the inverse + functions. + */ if (P < 0.05) { @@ -57,7 +63,7 @@ gsl_cdf_gamma_Pinv (const double P, const double a, const double b) else { double xg = gsl_cdf_ugaussian_Pinv (P); - double x0 = (xg < -sqrt (a)) ? a : sqrt (a) * xg + a; + double x0 = (xg < -0.5*sqrt (a)) ? a : sqrt (a) * xg + a; x = x0; } @@ -85,7 +91,7 @@ gsl_cdf_gamma_Pinv (const double P, const double a, const double b) double step1 = -((a - 1) / x - 1) * lambda * lambda / 4.0; double step = step0; - if (fabs (step1) < fabs (step0)) + if (fabs (step1) < 0.5 * fabs (step0)) step += step1; if (x + step > 0) @@ -140,7 +146,7 @@ gsl_cdf_gamma_Qinv (const double Q, const double a, const double b) else { double xg = gsl_cdf_ugaussian_Qinv (Q); - double x0 = (xg < -sqrt (a)) ? a : sqrt (a) * xg + a; + double x0 = (xg < -0.5*sqrt (a)) ? a : sqrt (a) * xg + a; x = x0; } @@ -168,7 +174,7 @@ gsl_cdf_gamma_Qinv (const double Q, const double a, const double b) double step1 = -((a - 1) / x - 1) * lambda * lambda / 4.0; double step = step0; - if (fabs (step1) < fabs (step0)) + if (fabs (step1) < 0.5 * fabs (step0)) step += step1; if (x + step > 0) _______________________________________________ Bug-gsl mailing list [email protected] http://lists.gnu.org/mailman/listinfo/bug-gsl
