Re: [Gretl-users] PDF of a multivariate normal
Thanks for your response, Jack.
Hehe. We're convergiung to the truth: here's the correct one
--- I think!
<hansl>
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 = dS * (2*$pi)^k
return exp(-0.5*sumr(U)) / sqrt(k)
end function
</hansl>
Ahh, ok, I see the issue. It works now, and also replicates gretl's
pdf-function in the univariate case as well.
>> 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.
Generating a sample of n independent Gaussian pseurdo-rvs with mean m
and covariance Sigma is rather trivial, in fact:
<hansl>
n = 20000
m = {1, 2, 3}
Sigma = {4,1,-1;1,4,1;-1,1,4}
A = mnormal(n, 3) * cholesky(Sigma)' .+ m
eval meanc(A)
eval mcov(A)
</hansl>
Artur