[Top][All Lists]

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
## [Bug-gsl] alternative implementation of gsl_ran_negative_binomial_pdf

**From**: |
Curtis Hash |

**Subject**: |
[Bug-gsl] alternative implementation of gsl_ran_negative_binomial_pdf |

**Date**: |
Thu, 02 Jun 2011 17:57:45 -0600 |

**User-agent**: |
Mozilla/5.0 (X11; U; Linux i686; en-US; rv:1.9.2.17) Gecko/20110424 Thunderbird/3.1.10 |

`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

**
**`nbinomial.c.patch`

*Description:* Text Data

**[Bug-gsl] alternative implementation of gsl_ran_negative_binomial_pdf**,
*Curtis Hash* **<=**