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

[Nils Bertschinger]

Hi Liam,

On Fri, 29 Jul 2011, Liam Healy wrote:

> 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.
Great, seems to work now. I did not pull the repository, but just checked 
your patch and applied it manually.

Best,

 	Nils

> 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>
>>
>