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]; |
| 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は要素が全て１の行列。これを使う事を思いつくのに |
| 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 |
/*苦労した割には、乱数でやった方が正確なので、あんまり報われない*/