help-gsl
[Top][All Lists]
Advanced

[Date Prev][Date Next][Thread Prev][Thread Next][Date Index][Thread Index]

[Help-gsl] Re: question about dirichlet lnpdf function (potential bug)


From: Rodney Sparapani
Subject: [Help-gsl] Re: question about dirichlet lnpdf function (potential bug)
Date: Mon, 13 Jul 2009 11:23:32 -0500
User-agent: Thunderbird 2.0.0.22 (X11/20090612)

per freem wrote:
 hi all,

i have been using the function for the pdf of the dirichlet distribution in
GSL (through a python interface.)

the function works correctly on most values, but appears to return nan when
it shouldn't on more extreme cases. for example, the pdf evaluated on the
set of values [0, 0, 1] with parameter settings [1/3, 1/3, 1/3] returns NaN,
even though [0, 0, 1] is a perfectly fine value for the dirichlet pdf that
should have non-zero probability. this is true for both the dirichlet_pdf
and dirichlet_lnpdf (log of pdf) functions.

in python notation, this is:

dirichlet_pdf([0, 0, 1], [1/3., 1/3., 1/3.])  (evaluates to NaN)
dirichlet_lnpdf([0, 0, 1], [1/3., 1/3., 1/3.])  (evaluates to NaN as well)

any idea why this is or how i can fix it?

thank you.

I decided to take a look:

#include <gsl/gsl_randist.h>

int main() {

  const double alpha[3]={0., 0., 0.}, theta[3]={1/3., 1/3., 1/3.};
  double x;
  x=gsl_ran_dirichlet_pdf(3, alpha, theta);

  return 0;
}

But, when I run this program, I get:

godzilla ~ % cc -g -L/usr/local/lib -I/usr/local/include \
-o dirichlet.out dirichlet.c -lgsl -lgslcblas -lm
godzilla ~ % dirichlet.out
gsl: gamma.c:1138: ERROR: domain error
Default GSL error handler invoked.
godzilla ~ % dbx dirichlet.out core
For information about new features see `help changes'
To remove this message, put `dbxenv suppress_startup_message 7.7' in your .dbxrc
Reading dirichlet.out
core file header read successfully
Reading ld.so.1
Reading libgsl.so.0.13.0
Reading libgslcblas.so.0.0.0
Reading libm.so.2
Reading libc.so.1
program terminated by signal ABRT (Abort)
0xfeaea067: __lwp_kill+0x0007: jae __lwp_kill+0x15 [ 0xfeaea075, .+0xe ]
Current function is gsl_error (optimized)
   47     abort ();
(dbx) where
  [1] __lwp_kill(0x1, 0x6), at 0xfeaea067
  [2] _thr_kill(0x1, 0x6), at 0xfeae57c4
  [3] raise(0x6), at 0xfea91da3
[4] abort(0x2, 0xfef29660, 0x8046e10, 0xfee1d178, 0xfef0ef2c, 0xfef0ec78), at 0xfea71bed =>[5] gsl_error(reason = ???, file = ???, line = ???, gsl_errno = ???) (optimized), at 0xfec7bdd5 (line ~47) in "error.c" [6] gsl_sf_lngamma_e(x = ???, result = ???) (optimized), at 0xfee1d178 (line ~1138) in "gamma.c" [7] gsl_sf_lngamma(x = ???) (optimized), at 0xfee1fdf4 (line ~1654) in "gamma.c" [8] gsl_ran_dirichlet_lnpdf(K = ???, alpha = ???, theta = ???) (optimized), at 0xfed9573b (line ~155) in "dirichlet.c" [9] gsl_ran_dirichlet_pdf(K = ???, alpha = ???, theta = ???) (optimized), at 0xfed954dc (line ~133) in "dirichlet.c"
  [10] main(), line 8 in "dirichlet.c"

It looks like we have log-gamma-function(0) which is 0.  However,
it is throwing an error.  Should it not just return 0?

Rodney





reply via email to

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