We believe there may be a more numerically stable way to write
gsl_ran_negative_binomial_pdf().
Test code:
#include <stdio.h>
#include <gsl/gsl_randist.h>
int main(int argc, char **argv) {
unsigned int k = 1015;
double n = 10150.0;
double p = n / (n + k);
printf("%.15f\n", gsl_ran_negative_binomial_pdf(k, p, n));
return 0;
}
The output of this test program with gsl-1.13+ is -nan.
However, using R:
dnbinom(1015, size=10150.0, p=(10150.0/(10150.0+1015.0)))
[1] 0.01193836
We believe this is because of the way P is calculated in
randist/nbinomial.c:
P = exp(f-a-b) * pow (p, n) * pow (1 - p, (double)k);
This method has problems with large values of n. Specifically,
exp(f-a-b) returns Inf.
Another way to express P is:
P = exp(f - a - b + n * log(p) + k * log(1 - p));
This will _always_ result in a smaller quantity passed to exp().
Attached is a patch against gsl-1.15/randist/nbinomial.c that implements
the suggested change.
The output of the test program with the patch is 0.011938361395305,
which is in agreement with R's implementation.
Josh Neil & Curt Hash
--- gsl-1.15-vanilla/randist/nbinomial.c 2010-12-26 10:57:08.000000000 -0700
+++ gsl-1.15/randist/nbinomial.c 2011-06-02 13:16:33.871584330 -0600
@@ -48,7 +48,7 @@ gsl_ran_negative_binomial_pdf (const uns
double a = gsl_sf_lngamma (n) ;
double b = gsl_sf_lngamma (k + 1.0) ;
- P = exp(f-a-b) * pow (p, n) * pow (1 - p, (double)k);
-
+ P = exp(f - a - b + n * log(p) + k * log(1 - p));
+
return P;
}
_______________________________________________
Bug-gsl mailing list
[email protected]
https://lists.gnu.org/mailman/listinfo/bug-gsl
