[Top][All Lists]
[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Bug-gsl] gsl_cdf_chisq_Pinv bug/feature
From: |
Bret Beck |
Subject: |
[Bug-gsl] gsl_cdf_chisq_Pinv bug/feature |
Date: |
Fri, 18 Apr 2008 10:24:57 -0700 |
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
The following code caused the problem.
-----------------------------------------------
---------------- Begin code --------------
-----------------------------------------------
#include <stdio.h>
#include <stdlib.h>
#include <gsl_rng.h>
#include <gsl_randist.h>
#include <gsl_cdf.h>
int main( int avgc, char **argv ) {
double nu = 1.5, P, v;
P = 1.93238e-01;
v = gsl_cdf_chisq_Pinv( P, nu );
printf( "P(%.17e) = %.17e\n", P, v );
P = 1.932383e-01;
v = gsl_cdf_chisq_Pinv( P, nu );
printf( "P(%.17e) = %.17e\n", P, v );
P = 1.93238145206123590e-01;
v = gsl_cdf_chisq_Pinv( P, nu);
printf( "P(%.17e) = %.17e\n", P, v );
exit( EXIT_SUCCESS );
}
-----------------------------------------------
----------------- End code ---------------
-----------------------------------------------
I looked into this and discovered the issue is in the routine
cdf/gammainv.c. When nu = 1.93238145206123590e-01 and nu = 0.75, the
number of allowed iterations at line 78 is too small. That is, I
change the line from
if (dP == 0.0 || n++ > 32)
to
if (dP == 0.0 || n++ > 40)
and it no longer fails. In fact, the problem is most likely in line
60, which is
double x0 = (xg < -sqrt (a)) ? a : sqrt (a) * xg + a;
For nu = 1.5 and around P = 1.93238145206123590e-01, xg is very close
to -sqrt( a ), and the guess x0 can jump from a = 0.75 to 0.
Thanks,
Bret Beck
- [Bug-gsl] gsl_cdf_chisq_Pinv bug/feature,
Bret Beck <=