[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
Re: [Bug-gsl] gsl_cdf_chisq_Pinv bug/feature
From: |
Brian Gough |
Subject: |
Re: [Bug-gsl] gsl_cdf_chisq_Pinv bug/feature |
Date: |
Tue, 29 Apr 2008 14:43:47 +0100 |
User-agent: |
Wanderlust/2.14.0 (Africa) Emacs/22.1 Mule/5.0 (SAKAKI) |
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)