function mrandomn, seed, covar, nrand, STATUS = status ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; ;+ ; NAME: ; MRANDOMN ; PURPOSE: ; Function to draw NRAND random deviates from a multivariate normal ; distribution with zero mean and covariance matrix COVAR. ; ; AUTHOR : Brandon C. Kelly, Steward Obs., Sept. 2004 ; ; INPUTS : ; ; SEED - The random number generator seed, the default is IDL's ; default in RANDOMN() ; COVAR - The covariance matrix of the multivariate normal ; distribution. ; OPTIONAL INPUTS : ; ; NRAND - The number of randomn deviates to draw. The default is ; one. ; OUTPUT : ; ; The random deviates, an [NRAND, NP] array where NP is the ; dimension of the covariance matrix, i.e., the number of ; parameters. ; ; OPTIONAL OUTPUT: ; STATUS - status of the Cholesky decomposition. If STATUS = 0 then ; the computation was successful. If STATUS > 0 then the ; input covariance matrix is not positive definite (see LA_CHOLDC), ; and MRANDOMN ; Note that if a STATUS keyword is supplied then no error message ; will be printed. ; REVISION HISTORY: ; Oct. 2013 -- Use LA_CHOLDC instead of CHOLDC to enable use of STATUS ; keyword. W. Landsman ;- ;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;;; if n_params() lt 2 then begin print, 'Syntax- Result = mrandomn( seed, covar, [nrand] , STATUS = )' return, 0 endif printerr = ~arg_present(errmsg) errmsg = '' ;check inputs and set up defaults if n_elements(nrand) eq 0 then nrand = 1 if size(covar, /n_dim) ne 2 then begin print, 'COVAR must be a matrix.' return, 0 endif np = (size(covar))[1] if (size(covar))[2] ne np then begin print, 'COVAR must be a square matrix.' return, 0 endif epsilon = randomn(seed, nrand, np) ;standard normal random deviates (NP x NRAND matrix) A = covar ;store covariance into dummy variable for input into TRIRED la_choldc, A, /double, status=status ;do Cholesky decomposition if status NE 0 then begin message,'Array is not positive definite, STATUS = ' + strtrim(status,2),/CON return,-1 endif for i = 0, np - 2 do A[i+1:*,i] = 0d ;Zero out upper triangular portion ;transform standard normal deviates so they have covariance matrix COVAR epsilon = A ## epsilon return, epsilon end