Re: [Bug-gsl] gsl_cdf_chisq_Pinv bug/feature

[Brian Gough]

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 <address@hidden>
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  <address@hidden>
+
+       * gammainv.c (gsl_cdf_gamma_Pinv, gsl_cdf_gamma_Qinv): restrict
+       the range of the gaussian approximation
+
 2008-02-20  Brian Gough  <address@hidden>
 
        * 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)