Wishart乱数の発生関数例

#include <oxstd.h>
#include <oxprob.h>
/////////////////////////////////////////////////
// Generation of Wishart Distribution using
// Bartlett's (1933) decompositon
// See Johnson, M.E. (1986)
// "Multivariate Statistical Simulation"
// Wiley. (page 204)
/////////////////////////////////////////////////
fMyWishart(const n,const mS){
// X~W(n,S), X:(pxp) symmetrx matrix
// f(X|n,S)
// = const *|X|^{(n-p-1)/2}*exp-0.5*tr(S^{-1}X)
// n: degrees of freedom
// S: (p x p) symmetric parameter matrix
// E(X_{ij})=n*s_{ij}, s_{ij}: (i,j) element of S.
decl ci,cj,cp,mA,mL,mT,mX;
cp=rows(mS);mT=zeros(cp,cp);
if(n<cp){print("Error! d.f. too small.");}
// generate A~W(n,I)
for(ci=0;ci<cp;++ci){
mT[ci][ci]=sqrt(ranchi(1,1,n-ci));
for(cj=0;cj<ci;++cj){mT[ci][cj]=rann(1,1);}
}
mA=mT*mT'; // A=TT' ~W(n,I)
mL=choleski(mS); // S=LL'(Choleski decomposition)
mX=mL*mA*mL'; // X=LAL'~W(n,S)
return(mX);
}

※あるいはVersion 3.3以降ではW(n,I)を発生させるranwishart(n,p)があるので、これを使ってmain文のなかで
mL=choleski(mS);mX=mL*ranwishart(n,p)*mL';
としてもよい(nは自由度(整数),pはXの次元(整数))。