モンテカルロシュミレーションにより円周率(パイ)を算出するプログラム。
”まぐれ”という不確実性科学の権威が書いた本を読んでいて簡単に
円周率を出せるというのでやってみた。
適当な円に密接する四角形を描く。ダーツゲームのようにその四角内に矢をあてる。
円の内部に入った矢と四角形内の矢すべての比率が円周率となる。
なんでかはわかりそうだったが、考えてたら頭がいたくなった・・・
とりあえず検証ということでrubyのTKでリアルタイムに描写させようとしたらはまった。
乱数発生の精度が低かった模様。最初は0〜100で乱数生成したので。
その後、精度を0〜10000にしたら良い数値が得られた。
3.14に近づくには100万回くらい待たなければだめみたい。
require "tkclass"
require "tk"
OFFSET=10
$oval=0
$sque=0
$disp_hankei=100
$hankei=10000
module MonteCarlo
def MonteCarlo.plot(max_x=$hankei*2, max_y=$hankei*2)
x = rand(max_x+1)
y = rand(max_y+1)
[x, y]
end
end
label = TkLabel.new do
text("0")
pack()
end
$c = TkCanvas.new('width'=>($disp_hankei+OFFSET)*2, 'height'=>($disp_hankei+OFFSET)*2).pack
Oval.new($c,OFFSET+0,OFFSET+0,OFFSET+$disp_hankei*2,OFFSET+$disp_hankei*2)
Rectangle.new($c,OFFSET+0,OFFSET+0,OFFSET+$disp_hankei*2,OFFSET+$disp_hankei*2)
def plot
x2 = 2
bx, by = MonteCarlo.plot($hankei*2, $hankei*2)
x, y = bx, by
x = bx - $hankei
y = by - $hankei
bx = bx / 100
by = by / 100
if (x*x + y*y) < $hankei*$hankei
$oval += 1
dot = Oval.new($c,OFFSET+bx,OFFSET+by,OFFSET+bx, OFFSET+by, "outline" => "red", "fill" => "red")
else
dot = Oval.new($c,OFFSET+bx,OFFSET+by,OFFSET+bx, OFFSET+by)
end
$sque += 1
end
count=0
TkAfter.new(1, -1, proc do
count+=1
plot
pai = $oval.to_f/$sque.to_f * 4.0
label.text("#{count.to_s} #{pai.to_s}")
end).start
Tk.mainloop
リアルタイムで描写しないのであれば乱数を少数で生成して円でなく扇がたで考え
シュミレーションするとあっというまによい精度でできた。
$oval = 0; $sque = 1
def plot
x, y = rand, rand
$oval += 1 if (x*x + y*y) < 1.0
$sque += 1
end
1000000.times {
plot
}
pai = $oval.to_f/$sque.to_f * 4.0
puts "pai #{pai.to_s}"
出力
pai 3.14324485675514
正直この本は難解な箇所が多々ありまが、理解できるとおもろいです。
