Re: [Gretl-users] PDF of a multivariate normal

matrix x = seq(-3,3.2,0.2)

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(a)univpm.it http://www2.econ.univpm.it/servizi/hpp/lucchetti ------------------------------------------------------- _______________________________________________ Gretl-users mailing list Gretl-users(a)lists.wfu.edu http://lists.wfu.edu/mailman/listinfo/gretl-users