[Bug-gsl] alternative implementation of gsl_ran_negative_binomial_pdf

[Curtis Hash]

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