MaTXでモンテカルロシミュレーション


MaTXは行列演算パッケージなので、基本的に演算はマトリクス変数でやるべきです。

MaTXはC言語っぽいスクリプトも書けますから

x=Z(1,5);
for (i=1;i<=5;i++){
  x(i)=rand();
}

なんていう様にforループを使いたいところですが、それはぐっと我慢して

x=rand(1,5)

と書きましょう。

forループを使ったら負け

とディスプレイの横にPostItで書いて張っておきましょう。(笑)
(実際、この手のMATLAB系のツールでは、forループは行列演算に比べるとかなり遅いです。)

MaTX Interpreter (matx)
Windows95/NT(Visual C++) version 5.0.5
last modified Fri Jan 28 18:28:19 JST 2000
Copyright (C) 1989-2000, Masanobu Koga

Send bugs and comments to matx@matx.org
Type 'quit' to exit, 'help' for functions, 'demo' for demonstration.


(以下、MaTXのプロンプトは省略します)

rand("seed"); //乱数初期化
x=rand(1,nmax); //nmax個の乱数ベクトル
y=rand(1,nmax); //nmax個の乱数ベクトル
r=x .* x + y .* y; // 要素毎の乗算。
   
t=(r .<=1.0 ); // 要素毎の論理演算 各要素<=1.0なら 1(真)そうでないなら 0(偽)
n=sum(t); // tの各要素につき t(i,j)==1 の個数を数える(tの要素を全て足す)
p=(n+0.0)/(nn * nn +0.0 )*4.0; //pの値を出力。 printf("PI=%f\n",p); の方がいいかも
p  
  p = 3.1404  

となって、π=3.1404と、結構イイ感じの値が。

 

/*******************************************************************/
乱数じゃなく、x,yを100等分して求めてみる。組み合わせで 100×100=10000個となる。
これ、forループでやれば、なんてことないんだが
forループを使ったら負け 
なので、
行列演算でやったら なんか結構複雑になっちまった。
/*******************************************************************/

nmax=100;  
x=linspace(0.0,1.0,nmax); //0〜1を nmax個で等分割したArray
y=linspace(0.0,1.0,nmax); //0〜1を等分割したArray
nn = length(x) //Arrayのサイズを念の為に求めなおす。(nmaxと同じ筈だけどね)
x2 = Matrix(x .* x); // x^2 を求めて、行列演算出来る様に、Matrix型にする。
y2= Matrix(y .* y)'; //こっちは ' をつけて列ベクトル化する
xx = ONE(nn,1) * x2; // ONEは要素が全て1の行列。これを使う事を思いつくのに
yy = y2 * ONE(1,nn); // 結構時間がかかった。(これで、行、又は列を全て同じにする)
r= xx + yy; // やっと x^2 + y^2 を求める。
   
t=(r .<= 1.0); // あとは一緒。
n=length(find(t)); // sum(t)でいいと思うんだが、本当はこうするべきなんだろな。
p=(n+0.0)/(nn * nn +0.0 )*4.0;  
p  
  p = 3.1156  

/*苦労した割には、乱数でやった方が正確なので、あんまり報われない*/