[GSLL-devel] Bug in dirichlet-log-pdf?

[Liam Healy]

Nils,

Thanks for the report.  I just pushed a fix to the git repository.  If you
can pull that and test it, please report if the results look right or not.
Otherwise Zach is planning a QL update in the next few days, so it should
show up then.

Liam

On Fri, Jul 29, 2011 at 8:49 AM, Nils Bertschinger <bertschi at mis.mpg.de>wrote:

> Hi,
>
> it appears that dirichlet-log-pdf is broken:
>
> CL-USER> (gsll:dirichlet-log-pdf #m(1.0d0 2.0d0) #m(0.3d0 0.7d0))
> 0.0d0
> CL-USER> (gsll:dirichlet-log-pdf #m(2.0d0 1.0d0) #m(0.7d0 0.3d0))
> -0.35667494393873245d0
>
> Obviously, the Dirichlet distribution is permutation invariant, so both
> results should be the same! By the way, both values are wrong. A straight
> forward implementation of the formula for the Dirichlet distribution gives:
>
> (defun dirichlet-log-pdf (alphas ps)
>  (let ((norm (loop for a across alphas
>                 sum (gsll:log-gamma a) into num
>                 sum a into denom
>                 finally (return (- num (gsll:log-gamma denom))))))
>    (- (reduce #'+ (map 'vector #'(lambda (p a) (* (- a 1.0d0) (log p)))
>                        ps alphas))
>       norm)))
>
> CL-USER> (dirichlet-log-pdf #(1.0d0 2.0d0) #(0.3d0 0.7d0))
> 0.33647223662121384d0
> CL-USER> (dirichlet-log-pdf #(2.0d0 1.0d0) #(0.7d0 0.3d0))
> 0.33647223662121384d0
> CL-USER>
>
> as it should be.
>
> Maybe I'm doing something wrong, when calling the gsl library function.
> Otherwise, I suspect there is a problem with passing arrays to the GSL.
> Broken memory layout?
>
> I'm using the following version of SBCL:
>
> CL-USER> (lisp-implementation-version)
> "1.0.45.0.debian"
> CL-USER> (software-version)
> "2.6.35-28-server"
> CL-USER> gsl:*gsl-version*
> "1.14"
>
> and my GSLL is from gsll-20101207-git as packaged in quicklisp.
>
> Any ideas of what is going on here?
>
> Thanks,
>
>        Nils
>
> ______________________________**_________________
> GSLL-devel mailing list
> GSLL-devel at common-lisp.net
> http://lists.common-lisp.net/**cgi-bin/mailman/listinfo/gsll-**devel<http://lists.common-lisp.net/cgi-bin/mailman/listinfo/gsll-devel>
>
-------------- next part --------------
An HTML attachment was scrubbed...
URL: <https://mailman.common-lisp.net/pipermail/gsll-devel/attachments/20110729/249248ee/attachment.html>