SICPを読む(92) 3.1.2 (2) - モンテカルロ法(1)
かなりゆっくり進んでおります。前回、M系列乱数を作ったので、モンテカルロ法を試しながら性能実験です。
SICPのモンテカルロ法がよくわからないので、とりあえず単純な奴を。
単純なモンテカルロ法を使って、円周率を量る
三平方の定理から、単位円上の円は
x^2 + y^2 = 1
となる。つまり
x^2 + y^2 <= 1
ならば、実数x,yは円の中にあるよね。
そこで、|x| <= 1, |y| <= 1の正方形中に、「でたらめに」点を打ってやると、
円の中に打たれた点 : 正方形の中に打たれた点 = 円の面積 : 正方形の面積
となり、面積の比が出る。
円の面積は πr^2で、単位円の半径はr = 1なので、単位円の面積=円周率。正方形の面積は4。つまり、
π/4 = 円の中に打たれた点/正方形の中に打たれた点
π = 4 * 円の中に打たれた点/正方形の中に打たれた点
となる。
スゲー単純な理論だけど、これで円周率が出るはずだ。
(define (monte-carilo/4 n) (letrec ((square (lambda (x) (* x x))) (iter (lambda (m) (if (zero? m) 0 (+ (if (<= (+ (square (rand)) (square (rand))) 1) 1 0) (iter (- m 1))))))) (* 4.0 (/ (iter n) n))))
テスト。
(monte-carilo/4 10) ; 3.2 (monte-carilo/4 100) ; 3.08 (monte-carilo/4 1000) ; 3.168 (monte-carilo/4 10000) ; 3.1256 (monte-carilo/4 100000) ; 3.13932
・・・。
「全く収束しない」
1万回くらいやると3.1が決定できるくらいなので、3.141位まで出そうと思うだけで泣けますね。
