Multivariate normal pdf code?

[ignacio82]

I'm looking for Fortran code that calculates the multivariate normal pdf.

If someone is willing to share or point me in the right direction I would appreciate it.

Thanks!

[ignacio82]

I have code running that do this (see attachment).

But, if someone has something that is faster it would be great.

Thanks again!

[jimdempseyatthecove]

mecej4,

This is a great example of why FORTRAN is still used.

rtdet_sigma=product([(sigma(i,i),i=1,k)])
mv_normal_pdf= max(exp(inside_exp)/(c*rtdet_sigma),1d-16)

Are somewhat hard to do in C++

Jim Dempsey

[mecej4]

ignacio82 wrote:
But, if someone has something that is faster it would be great.

Here is my offering, which uses MKL to compute the Cholesky factorization of the covariance matrix for the "nondegenerate" case where the matrix is positive-definite. CAUTION: please test with at least a few example cases before relying on this lightly tested code. Note also that the calls to MKL cause the input argument Σ to be overwritten by the factor, so if you want Σ to be preserved, you should declare a local array, copy Σ to it, and use the copy of Σ in the call.

All that I have done to your code is to replace several lines by calls to LAPACK and BLAS routines. As the forum software mangles the source code when displaying it inline, especially so the indentation, I have attached the file, too.

[fortran]

module stats
implicit none
contains
double precision function mv_normal_pdf(k,x,x_mean,sigma)
use mkl95_lapack, only : potrf,trtrs
use mkl95_blas, only : dot
implicit none
integer, intent(in) :: k
double precision, intent(in) :: x(k), x_mean(k)
double precision, intent(in out) :: sigma(k,k)
double precision :: rtdet_sigma, xminusmean(k)
doubleprecision :: c, inside_exp
integer :: i,info

call potrf( sigma, 'L' ,info )
if(info /= 0)then
write(*,*)' INFO from POTRF = ',info
stop
endif
rtdet_sigma=product([(sigma(i,i),i=1,k)])
xminusmean=x-x_mean
call trtrs(sigma,xminusmean,'L','N','N',info)
if(info /= 0)then
write(*,*)' INFO from TRTRS = ',info
stop
endif
inside_exp=-0.5d0*dot(xminusmean,xminusmean)
c = (atan(1d0)*8d0)**(k/2d0)
mv_normal_pdf= max(exp(inside_exp)/(c*rtdet_sigma),1d-16)
end function mv_normal_pdf

end module stats

[/fortran]