モンテカルロ法のアニメーション
乱数を使った数値処理は色々ある。暗号の世界では予測出来ないような乱数が重要だったような？
また、シェル芸本に数値データを扱う記事があり、そこでモンテカルロ法で円周率を求める方法が記載されていた。 gihyo.jp
面積を求める時に乱数座標で打ち抜くモンテカルロ法ってのがある。モンテカルロ法では、一様乱数の質が良いことが求められるはず？じゃあ実際にプロットされる乱数座標がどんな感じなのか、アニメーションGIFを作ってみることにした。
乱数というかランダムビット生成には/dev/urandomという乱数生成仮想デバイスを使うのがお手軽だが、LinuxやMacでは読み出し速度があんまり速くない。実験した限りではFreeBSDが比較的早い。
awkにrand()という0から1の範囲で乱数を生成する関数がある。/dev/urandomから読み出すより処理が速かったので、今回はawkのrand()関数を使うことにする。
ということで、100万の擬似一様乱数座標をプロットし、円周率が収束していくアニメーション作成のシェルスクリプトを作ってみた。実行環境は下記の通り。
- OSX El Capitan
- GNU coreutils
- GNU awk
- ImageMagick
- gnuplot
- Open usp Tukubai
作成したシェルスクリプトは下記の通り。実行するとanimate.gifという名前のアニメーションGIFが出来る。gnuplotスクリプトを、Tukubaiのmojihameコマンドを使って作成している。
pi_gif.bash
#!/bin/bash gseq 1000000 | gawk 'BEGIN{srand()}{x=rand()-0.5;y=rand()-0.5;if(sqrt(x^2+y^2)<0.5)n++; print NR,x,y,4*n/NR}' > data.txt cat << EOF > template.gp set terminal png size 800,800 set nokey set xr[0:0.5] set yr[0:0.5] set size square set parametric set samples 1000 set title font 'Helvetica,20' LABEL set output "%1.png" set title "N=%1 pi=%2"; plot 'data.txt' every ::0::%3 using 2:3 with points pointsize 0.2, [0:2*pi] cos(t)/2,sin(t)/2 linewidth 0.5; LABEL EOF sed -n '10p;100p;1000p;10000p;100000p;$p' data.txt | awk '{print $1,$NF,$1-1}' | mojihame -lLABEL template.gp - > plot.gp gnuplot plot.gp convert -delay 300 *.png animate.gif
$ ./pi_gif.bash $ ls -l *.png *.gif -rw-r--r--@ 1 tashiro staff 5406 12 20 17:43 10.png -rw-r--r--@ 1 tashiro staff 5890 12 20 17:43 100.png -rw-r--r--@ 1 tashiro staff 7988 12 20 17:43 1000.png -rw-r--r--@ 1 tashiro staff 20838 12 20 17:43 10000.png -rw-r--r--+ 1 tashiro staff 67217 12 20 17:43 100000.png -rw-r--r--+ 1 tashiro staff 44507 12 20 17:43 1000000.png -rw-r--r--@ 1 tashiro staff 136680 12 20 17:43 animate.gif
乱数座標の数と円周率を、画像上部のタイトルに表示してある。こんな感じのアニメーションになった。
乱数座標は密度のバラツキが多少ある。もっと綺麗にバラついた一様乱数座標があると、精度の高い円周率を求めることが出来るのかもしれない。
円周率は (円の中に入った座標数/全体の座標数)*4 で計算される。
