bug-gsl
[Top][All Lists]
Advanced

[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)))
 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

reply via email to

 [Prev in Thread] Current Thread [Next in Thread]