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


モンテカルロっていうのは、乱数を多数発生させて、 …と書いていて気がついたのですが、SCILABって乱数のSeed(種)を例えば時間を使って 乱数にするっていうのはどうすりゃいいんだか?

rand("seed",n) でnにできるんだが、これだと、毎回同じ乱数列しか得られないので 例えば普通は
rand("seed",second) とかして秒なり分をSeedにするんだが、どういうわけか gettime()では年と月しか出ないし。

-->getdate()
ans =
! 2000. 10. 0. 0. 0. 0. 0. 0. 0. !

こりゃ困ったね(笑)。まぁいいか。


SCILABでモンテカルロシミュレーションで円周率を求めよう。

(どうやって円周率を求めるか?の前振りはここ)

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

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

for i=1:5 , x(i)=rand(); , end

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

x=rand(1,5)

と書きましょう。

forループを使ったら負け

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


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

nmax=10000; //nmax組の乱数で行う。
x=rand(1,nmax); //nmax個の乱数ベクトル
y=rand(1,nmax); //nmax個の乱数ベクトル
r=x .* x + y .* y; // 要素毎の乗算。 x * x じゃなくて、 x * xであることに注意
   
t = find( r<=1.0 );

// 各要素<=1.0である要素の番号をリストアップする

例:

-->A=[1 2 3 4 5 4 3 2 1];
-->find( A<=3 )
ans =
! 1. 2. 3. 7. 8. 9. !

n= length(t); // tの要素数を数える。
p=(n)/(nmax)*4.0; //pの値を出力。 printf("PI=%f\n",p); の方がいいかも
p  
  p = 3.1564  

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

 

/*******************************************************************/
乱数じゃなく、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 = (x .* x); // x^2 を求めて、行列演算出来る様に、Matrix型にする。
y2= (y .* y)'; //こっちは ' をつけて列ベクトル化する
xx = ones(nn,1) * x2; // onesは要素が全て1の行列。これを使う事を思いつくのに
yy = y2 * ones(1,nn); // 結構時間がかかった。(これで、行、又は列を全て同じにする)
r= xx + yy; // やっと x^2 + y^2 を求める。
   
t=find(r <= 1.0); // あとは一緒。
n=length(t);  
p=(n)/(nmax * nmax) *4.0;  
p  
  p = 3.1156  

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


return