On Tue, 16 Feb 2016, Artur T. wrote:

Ok, I changed the definition of the constant k. However, there is still
some difference to the outcome using MATLAB...

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>

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>

-------------------------------------------------------
  Riccardo (Jack) Lucchetti
  Dipartimento di Scienze Economiche e Sociali (DiSES)

  Università Politecnica delle Marche
  (formerly known as Università di Ancona)

  r.lucchetti@univpm.it
  http://www2.econ.univpm.it/servizi/hpp/lucchetti
-------------------------------------------------------

_______________________________________________
Gretl-users mailing list
Gretl-users@lists.wfu.edu
http://lists.wfu.edu/mailman/listinfo/gretl-users