Ox（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の次元(整数))。

#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