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

[Brian Gough]

At Sun, 12 Jul 2009 13:09:37 -0400,
per freem wrote:
> 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?

Can you confirm which way round you're using the parameters?  Based on
the description in your email, the order of the python call (theta,
alpha) seems the opposite to the order of the GSL arguments (alpha,
theta).

If the parameters are alpha=(1/3,1/3,1/3) then the pdf of (0,0,1)
should be infinite since it is proportional to theta[0]^-2/3
theta[1]^-2/3 theta[2]^-2/3 and this is what I get.

If it's the other way round there is a problem with alpha values of
zero, and the patch below might fix it.

-- 
Brian Gough
(GSL Maintainer)

Support freedom by joining the FSF 
http://www.fsf.org/associate/support_freedom/join_fsf?referrer=37

--- a/randist/dirichlet.c
+++ b/randist/dirichlet.c
@@ -156,7 +156,7 @@ gsl_ran_dirichlet_lnpdf (const size_t K,
 
   for (i = 0; i < K; i++)
     {
-      log_p -= gsl_sf_lngamma (alpha[i]);
+      log_p -= (alpha[i] == 0) ? GSL_POSINF : gsl_sf_lngamma (alpha[i]);
     }
 
   return log_p;