bugGNU Scientific Library - Bugs: bug #45924, Bug in the inverse beta function...

 
 

You are not allowed to post comments on this tracker with your current authentication level.

bug #45924: Bug in the inverse beta function gsl_cdf_beta_Pinv, and suggested fix

Submitter:  Patrick Alken <psa>
Submitted:  Fri 11 Sep 2015 10:06:25 AM UTC
   
 
Category:  Runtime error Severity:  3 - Normal
Operating System:  Status:  None
Assigned to:  None Open/Closed:  Open
Release: 

Fri 11 Sep 2015 10:06:25 AM UTC, original submission:  

from jorg =at= google =dot= com

For low or high values of P, the inverse beta function fails to converge.
Here are some inputs and results (the boost inverse-beta function is
included for context.)

a: 8000, b:2000, p: 0.995,   boost: 0.810189,  gsl:nan
a: 8000, b:2000, p: 0.005,   boost: 0.789585,  gsl:nan
a: 8000, b:2000, p: 0.99995, boost: 0.815275,  gsl:nan
a: 8000, b:2000, p: 5e-05,   boost: 0.78416,   gsl:nan
a: 630,  b:9370, p: 0.995,   boost: 0.0694215, gsl:nan
a: 630,  b:9370, p: 0.99995, boost: 0.0728637, gsl:nan
a: 5000, b:5000, p: 0.995,   boost: 0.512877,  gsl:nan
a: 5000, b:5000, p: 0.005,   boost: 0.487123,  gsl:nan
a: 5000, b:5000, p: 0.99995, boost: 0.519446,  gsl:nan
a: 5000, b:5000, p: 5e-05,   boost: 0.480554,  gsl:nan

As I looked into this bug, what appears to be going on is that the call to
bisect does an early return, thinking that the initial guess is good
enough.  But the early guess is too far away for the rest of the routine to
converge.

If I change these two lines:

  /* Do bisection to get closer */
  x = bisect (x, P, a, b, 0.01, 0.01);

To this:

  /* Do bisection to get closer */
  if (P < 0.1) {
    x = bisect (x, P, a, b, P/10, P/10);
  } else {
    x = bisect (x, P, a, b, 0.01, 0.01);
  }

Then the problem goes away.

Boilerplate follows:
Version number of GSL: 1.16, however the latest sources downloaded from git
appear to have the same problem.
Hardware and OS: HP Z620, Ubuntu Trusty
Compiler used: either gcc or clang, optimized or not.
Short program which exercises the bug:

(Code is in C++11)

#include <iostream>
#include "third_party/gsl/gsl_1_16/cdf/gsl_cdf.h"

int main(int argc, char *argv[]) {
  std::cout << "Hello World\n";
  struct Test {
    double a;
    double b;
    double p;
    double x;
  } test_cases[] = {
                 {8000, 2000, 0.995,   0.810189},
                 {8000, 2000, 0.005,   0.789585},
                 {8000, 2000, 0.99995, 0.815275},
                 {8000, 2000, 0.00005, 0.78416},
                 { 630, 9370, 0.995,   0.0694215},
                 { 630, 9370, 0.99995, 0.0728637},
                 {5000, 5000, 0.995,   0.512877},
                 {5000, 5000, 0.005,   0.487123},
                 {5000, 5000, 0.99995, 0.519446},
                 {5000, 5000, 0.00005, 0.480554},
               };
  // beta function
  for (auto t : test_cases) {
    std::cout << "a=" << t.a << " b=" << t.b << " x=" << t.x << " gsl=" <<
gsl_cdf_beta_P(t.x, t.a, t.b) << "\n";
  }
  // inverse beta function
  for (auto t : test_cases) {
    std::cout << "a=" << t.a << " b=" << t.b << " p=" << t.p << " gsl=" <<
gsl_cdf_beta_Pinv(t.p, t.a, t.b) << "\n";
  }
}

Patrick Alken <psa>
Group administrator

 

(Note: upload size limit is set to 16384 kB, after insertion of the required escape characters.)

No files currently attached

 

Depends on the following items: None found

Items that depend on this one: None found

 

Carbon-Copy List
  • -email is unavailable- added by psa (Submitted the item)
  •  

    No changes have been made to this item

    Back to the top

    Powered by Savane 3.13-758e.
    Corresponding source code