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.
I have code running that do this (see attachment).
But, if someone has something that is faster it would be great.
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.
double precision function mv_normal_pdf(k,x,x_mean,sigma)
use mkl95_lapack, only : potrf,trtrs
use mkl95_blas, only : dot
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
if(info /= 0)then
write(*,*)' INFO from TRTRS = ',info
c = (atan(1d0)*8d0)**(k/2d0)
end function mv_normal_pdf
end module stats
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++