Code: Select all
'command: wish(S, h, n) -- generates multivariate Wishart random variables
'Inputs: h --- m-by-m scale matrix
' n --- scalar degree of freedom
'Outputs: S = m-by-m matrix draw from the Wishart distribution
'Note: used Bartlett decomposition
'Example:
matrix(2,2) s
sym(2) h = @identity(2)
scalar n = 5
call wish(s,h,n)
subroutine local wish(matrix S, sym h, scalar n)
!p = @rows(h)
matrix T = @zeros(!p,!p)
for !i=1 to !p
for !j=1 to !i
if !i = !j then
!ntemp = @mrnd(1,1)
!chirnd = @qchisq(!ntemp,n-!i+1)
T(!i,!j) = @sqrt(!chirnd)
else
!normrnd = @mnrnd(1,1)
T(!i,!j) = !normrnd
endif
next
next
matrix A = T*@transpose(T)
matrix Y = @cholesky(h)
S = @transpose(Y)*A*Y
endsub