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

[Nils Bertschinger]

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