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

On Mon, 15 Feb 2016, Artur T. wrote: OK, you win, dammit! <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 = (2*$pi)^k / sqrt(dS) return k .* exp(-0.5*sumr(U)) end function nulldata 10000 series x = normal() series y = 2*x + normal() matrix A = {5, 1; 1, 1} matrix C = {y, x} matrix m = zeros(2,1) series f = mnormal_pdf(C, m, 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 -------------------------------------------------------