Hi,
The function gsl_ran_beta_pdf contains the following section,
which is causing incorrect overflows for x == 0 or 1 and when
a and b are large enough:
(at randist/beta.c:53)
...
double gab = gsl_sf_lngamma (a + b);
double ga = gsl_sf_lngamma (a);
double gb = gsl_sf_lngamma (b);
if (x == 0.0 || x == 1.0)
{
p = exp (gab - ga - gb) * pow (x, a - 1) * pow (1 - x, b - 1);
}
else
...
For example:
(gdb) call gsl_ran_beta_pdf(1,100,100)
$26 = 0
(gdb) call gsl_ran_beta_pdf(0,100,100)
$27 = 0
(gdb) call gsl_ran_beta_pdf(1,1000,1000)
$28 = -nan(0x8000000000000)
(gdb) call gsl_ran_beta_pdf(0,1000,1000)
$29 = -nan(0x8000000000000)
but the correct result is always 0 when x == 0 or 1 and a, b > 1.
The reason is that the call to exp() overflows when a and b are
large enough. However, it's unnecessary to call exp() when either
of the two calls to pow() will return 0.
I'm attaching a minimal patch which simply checks for a, b > 1
and explicitly returns 0 in that case.
Mike
--- beta.c 2010-03-10 02:57:13.000000000 -0800
+++ beta-proposed.c 2010-07-15 15:12:55.941011520 -0700
@@ -56,7 +56,14 @@
if (x == 0.0 || x == 1.0)
{
- p = exp (gab - ga - gb) * pow (x, a - 1) * pow (1 - x, b - 1);
+ if (a > 1.0 && b > 1.0)
+ {
+ p = 0.0;
+ }
+ else
+ {
+ p = exp (gab - ga - gb) * pow (x, a - 1) * pow (1 - x, b - 1);
+ }
}
else
{
_______________________________________________
Bug-gsl mailing list
[email protected]
http://lists.gnu.org/mailman/listinfo/bug-gsl
