top page > computer > haskell > web_lecture > for_programmer > montecarlo_coding.html
更新日:
文責: 重城良国

モンテカルロ法のコーディング

montecarlo.hs

正方形内のランダムな点

正方形内のランダムな点を生成する式は以下のようだった。

uncurry zip . (randomRs (-1, 1) *** randomRs (-1, 1)) . split $ mkStdGen 8

乱数の種の生成のための整数を引数としいろいろな系列のランダムな点の列を返す関数とする。

points :: Int -> [(Double, Double)]
points = uncurry zip . (randomRs (-1, 1) *** randomRs (-1, 1)) . split . mkStdGen

円内の点であることのチェック

原点を中心とする半径1の円のなかとは原点からの距離が1以下ということだ。

inCircle :: (Double, Double) -> Bool
inCircle (x, y) = x ^ 2 + y ^ 2 <= 1

距離の2乗が1以下と考えている。

円内の点のみを取り出す

全体の点の数を指定して円内の点のみをとりだす関数だ。

inCirclePoints :: Int -> Int -> [(Double, Double)]
inCirclePoints g n = filter inCircle . take n $ points g

円周率の推測値を出す

円内の点の数を全体の点の数で割ったものに正方形の面積である4をかける。

[正方形の面積] * [円内の点の数] / [正方形内の点の数]

ランダムな点の系列と点の数を指定し円周率の推測値を計算する。

guessPi :: Int -> Int -> Double
guessPi g n = 4 * fromIntegral (length $ inCirclePoints g n) / fromIntegral n

試してみる

% ghci montecarlo.hs
*Main> guessPi 8 100
2.96
*Main> guessPi 8 1000
3.096
*Main> guessPi 8 10000
3.134
*Main> guessPi 8 100000
3.1434
*Main> guessPi 8 1000000
3.143328
*Main> guessPi 9 1000000
3.142572

グラフ

横軸が点の数、縦軸が円周率の推測値だ。収束が遅いので横軸は対数だ。gの値により異なる系列のランダム列となる。

[円周率へ収束していくグラフ]

πに収束しているように見える。

まとめ

モンテカルロ法で円周率を求めた。関数型言語の魂である「再帰」をまだ学んでいないがリストだけでここまでできる。

「疑似乱数列の生成」へもどる 「構文: let ... in ...」へ

正当なCSSです! HTML5 Powered with CSS3 / styling, and Semantics