読者です 読者をやめる 読者になる 読者になる

日々之迷歩

世の中わからんことだらけ

ITが複雑で難しくなっていく様に翻弄される日々です

モンテカルロ法のアニメーション

乱数を使った数値処理は色々ある。暗号の世界では予測出来ないような乱数が重要だったような?

また、シェル芸本に数値データを扱う記事があり、そこでモンテカルロ法で円周率を求める方法が記載されていた。 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 で計算される。

f:id:papiro:20151220175021g:plain