[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]
[Bug-gsl] [bug #45924] Bug in the inverse beta function gsl_cdf_beta_Pin
From: |
Patrick Alken |
Subject: |
[Bug-gsl] [bug #45924] Bug in the inverse beta function gsl_cdf_beta_Pinv, and suggested fix |
Date: |
Fri, 11 Sep 2015 10:06:26 +0000 |
User-agent: |
Mozilla/5.0 (X11; Ubuntu; Linux x86_64; rv:40.0) Gecko/20100101 Firefox/40.0 |
URL:
<http://savannah.gnu.org/bugs/?45924>
Summary: Bug in the inverse beta function gsl_cdf_beta_Pinv,
and suggested fix
Project: GNU Scientific Library
Submitted by: psa
Submitted on: Fri 11 Sep 2015 10:06:25 AM GMT
Category: Runtime error
Severity: 3 - Normal
Operating System:
Status: None
Assigned to: None
Open/Closed: Open
Release:
Discussion Lock: Any
_______________________________________________________
Details:
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";
}
}
_______________________________________________________
Reply to this item at:
<http://savannah.gnu.org/bugs/?45924>
_______________________________________________
Message sent via/by Savannah
http://savannah.gnu.org/
[Prev in Thread] |
Current Thread |
[Next in Thread] |
- [Bug-gsl] [bug #45924] Bug in the inverse beta function gsl_cdf_beta_Pinv, and suggested fix,
Patrick Alken <=