/usr/share/gnudatalanguage/astrolib/mrandomn.pro is in gdl-astrolib 2018.02.16+dfsg-1.
This file is owned by root:root, with mode 0o644.
The actual contents of the file can be viewed below.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 | 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
|