Multivariate normal pdf code?

Multivariate normal pdf code?

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!

4 posts / 0 new
Last post
For more complete information about compiler optimizations, see our Optimization Notice.

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

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

Thanks again!

Attachments: 

AttachmentSize
Download stats-module.f907.72 KB

Citação:

ignacio82 escreveu:
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.

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

Attachments: 

AttachmentSize
Download stats-mkl.f901.28 KB

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

www.quickthreadprogramming.com

Leave a Comment

Please sign in to add a comment. Not a member? Join today