Nils,<br><br>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.<br>

<br>Liam<br><br><div class="gmail_quote">On Fri, Jul 29, 2011 at 8:49 AM, Nils Bertschinger <span dir="ltr"><<a href="mailto:bertschi@mis.mpg.de">bertschi@mis.mpg.de</a>></span> wrote:<br><blockquote class="gmail_quote" style="border-left: 1px solid rgb(204, 204, 204); margin: 0pt 0pt 0pt 0.8ex; padding-left: 1ex;">

Hi,<br>
<br>
it appears that dirichlet-log-pdf is broken:<br>
<br>
CL-USER> (gsll:dirichlet-log-pdf #m(1.0d0 2.0d0) #m(0.3d0 0.7d0))<br>
0.0d0<br>
CL-USER> (gsll:dirichlet-log-pdf #m(2.0d0 1.0d0) #m(0.7d0 0.3d0))<br>
-0.35667494393873245d0<br>
<br>
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:<br>


<br>
(defun dirichlet-log-pdf (alphas ps)<br>
  (let ((norm (loop for a across alphas<br>
                 sum (gsll:log-gamma a) into num<br>
                 sum a into denom<br>
                 finally (return (- num (gsll:log-gamma denom))))))<br>
    (- (reduce #'+ (map 'vector #'(lambda (p a) (* (- a 1.0d0) (log p)))<br>
                        ps alphas))<br>
       norm)))<br>
<br>
CL-USER> (dirichlet-log-pdf #(1.0d0 2.0d0) #(0.3d0 0.7d0))<br>
0.33647223662121384d0<br>
CL-USER> (dirichlet-log-pdf #(2.0d0 1.0d0) #(0.7d0 0.3d0))<br>
0.33647223662121384d0<br>
CL-USER><br>
<br>
as it should be.<br>
<br>
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?<br>
<br>
I'm using the following version of SBCL:<br>
<br>
CL-USER> (lisp-implementation-version)<br>
"1.0.45.0.debian"<br>
CL-USER> (software-version)<br>
"2.6.35-28-server"<br>
CL-USER> gsl:*gsl-version*<br>
"1.14"<br>
<br>
and my GSLL is from gsll-20101207-git as packaged in quicklisp.<br>
<br>
Any ideas of what is going on here?<br>
<br>
Thanks,<br>
<br>
        Nils<br>
<br>
______________________________<u></u>_________________<br>
GSLL-devel mailing list<br>
<a href="mailto:GSLL-devel@common-lisp.net" target="_blank">GSLL-devel@common-lisp.net</a><br>
<a href="http://lists.common-lisp.net/cgi-bin/mailman/listinfo/gsll-devel" target="_blank">http://lists.common-lisp.net/<u></u>cgi-bin/mailman/listinfo/gsll-<u></u>devel</a><br>
</blockquote></div><br>