Random draw from multivariate Normal
Posted: Wed Jul 02, 2014 10:44 am
Below is a subroutine to produce a sample from the specified multivariate normal distribution. Comments are welcomed.
Code: Select all
'command: vector(k) y; vector(k) mean; sym(k) cov; mvrnd(y, mean, cov)
will generate a multivariate normal random variate y
'Inputs: mean --- k-by-1 mean vector
' cov --- k-by-k symmetric covariance
'Outputs: y = k-by-1 normal draw from the multivariate normal distribution
subroutine local mvnrnd(vector y, vector mean, sym cov)
!n = @rows(cov)
if !n = 1 then
y = mean + @transpose(@sqrt(cov))*@mnrnd(!n,1)
else if !n > 1 then
matrix cc = @cholesky(cov)
y = mean + @transpose(cc)*@mnrnd(!n,1)
endif
endif
endsub