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

Attachment: nbinomial.c.patch
Description: Text Data

reply via email to

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