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

Clive, Artur's version of gretl is newer than yours or mine (though I've since updated). The integer requirements for seq() have been removed. PS Sent from my iPad On Feb 16, 2016, at 9:33 PM, Clive Nicholas <clivelists@googlemail.com<mailto:clivelists@googlemail.com>> wrote: Thanks Peter, that sorted it. Do remember that this wasn't my code, so how Artur's original code worked for him I've no clue. C On 16 February 2016 at 18:50, Summers, Peter <psummers@highpoint.edu<mailto:psummers@highpoint.edu>> wrote: Clive, The arguments to seq() all have to be integers. Try replacing that line with “matrix x = seq(-30,32,2)./10” PS From: gretl-users-bounces@lists.wfu.edu<mailto:gretl-users-bounces@lists.wfu.edu> [mailto:gretl-users-bounces@lists.wfu.edu<mailto:gretl-users-bounces@lists.wfu.edu>] On Behalf Of Clive Nicholas Sent: Tuesday, February 16, 2016 1:27 PM To: r.lucchetti@univpm.it<mailto:r.lucchetti@univpm.it>; Gretl list Subject: Re: [Gretl-users] PDF of a multivariate normal Jack, Artur's code chokes for me when I run your function code: gretl version 2016a-git Current session: 2016-02-16 18:17 Invalid argument Error executing script: halting I wonder why? My -gretl-'s fully up to date! C On 16 February 2016 at 15:26, Riccardo (Jack) Lucchetti <r.lucchetti@univpm.it<mailto:r.lucchetti@univpm.it>> wrote: 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<mailto:r.lucchetti@univpm.it> http://www2.econ.univpm.it/servizi/hpp/lucchetti ------------------------------------------------------- _______________________________________________ Gretl-users mailing list Gretl-users@lists.wfu.edu<mailto:Gretl-users@lists.wfu.edu> http://lists.wfu.edu/mailman/listinfo/gretl-users -- Clive Nicholas "My colleagues in the social sciences talk a great deal about methodology. I prefer to call it style." -- Freeman J. Dyson _______________________________________________ Gretl-users mailing list Gretl-users@lists.wfu.edu<mailto:Gretl-users@lists.wfu.edu> http://lists.wfu.edu/mailman/listinfo/gretl-users -- Clive Nicholas "My colleagues in the social sciences talk a great deal about methodology. I prefer to call it style." -- Freeman J. Dyson _______________________________________________ Gretl-users mailing list Gretl-users@lists.wfu.edu<mailto:Gretl-users@lists.wfu.edu> http://lists.wfu.edu/mailman/listinfo/gretl-users