#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の次元(整数))。
#include <oxstd.h>
fTest1(const x){
decl y;
y = x;println("(2) x=",y);
}
fTest2(const adX){
println("(4) x=",adX[0]);
adX[0]=1;
println("(5) x=",adX[0]);
}
main(){
decl x;
x = 0; println("(1) x=",x);
println("Function without using address");
fTest1(x); println("(3) x=",x);
//
println("Function using address");
fTest2(&x);
println("(6) x=",x);
}
出力結果
(1) x=0
Function without using address
(2) x=0
(3) x=0
Function using address
(4) x=0
(5) x=1
(6) x=1
#include <oxstd.h>
fSquare(const y){
decl y2;
y2 = y^2;
return y2;
}
main(){
decl x,x2;
x = 2;
x2 = fSquare(x);
println("x=",x,", x^2=",x2);
}
出力結果
x=2, x^2=4