Ok, I changed the definition of the constant k. However, there is still
some difference to the outcome using MATLAB... Attached is also a
comparison --in the univarite case -- using gretl's pdf-function.
Actually, I am not sure whether these differences are due to numerical
issues or whether there is still some issue in the "mnormal_pdf"-function
<hansl>
clear
set echo off
set messages off
function matrix mnormal_pdf(matrix x, matrix m, matrix Sigma)
# x should be n x k (n k-dimensional vectors)
# m should be a k-vector of means
# Sigma should be a kxk pd matrix
scalar k = cols(x)
scalar n = rows(x)
if k != rows(m)*cols(m)
return mshape(NA, n, 1)
endif
scalar dS = det(Sigma)
if dS <= 0
return mshape(NA, n, 1)
endif
matrix U = x .- vec(m)'
matrix U = U .* (U/Sigma)
#scalar k = (2*$pi)^k / sqrt(dS) # OLD
scalar k = 1/(sqrt(dS)*(2*$pi)^k)
return k * exp(-0.5*sumr(U))
end function
# Bivariate Example
matrix mu = zeros(2,1)
matrix Sigma = {.25, .3; .3, 1}
matrix x = seq(-3,3.2,0.2)
matrix X = x' ~ x'
matrix f = mnormal_pdf(X, mu, Sigma)
f = f ~ x'
gnuplot 1 2 --matrix=f --with-lines --output=display
# CHECK: If Sigma is a scalar the problem reduces to the univariate #normal
matrix mu = zeros(1,1)
matrix Sigma = {1}
matrix f2 = mnormal_pdf(X[,1], mu, Sigma)
matrix f3 = pdf(N,X[,1])
f = f2 ~ f3 ~ x'
gnuplot 1 2 3 --matrix=f --with-lines --output=display
</hansl>
Actually I was a little surprised about Jack's suggestion also
with
respect to what it does. I thought a function is wanted that also does
the RNG inside.
No, it is really just about computing the pdf.
Artur