モンテカルロでパイ

モンテカルロシュミレーションにより円周率(パイ)を算出するプログラム。
”まぐれ”という不確実性科学の権威が書いた本を読んでいて簡単に
円周率を出せるというのでやってみた。

適当な円に密接する四角形を描く。ダーツゲームのようにその四角内に矢をあてる。
円の内部に入った矢と四角形内の矢すべての比率が円周率となる。
なんでかはわかりそうだったが、考えてたら頭がいたくなった・・・

とりあえず検証ということで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

正直この本は難解な箇所が多々ありまが、理解できるとおもろいです。