疑似乱数に状態なんていらない

State モナドと疑似乱数で書いたように、遅延評価が利用できる言語では、無限数列が扱えるので、疑似乱数を使う際に状態を持たなくてもよい。その一例として、モンテカルロ法による円周率の近似を挙げてみる。

XY 平面に単位円を考える。

radius :: Double
radius = 1.0

この円がぴったり収まる大きさ1の正方形を描く。ここで、第一象限のみを考える。正方形のうち、第一象限にある部分の面積は、1/4。第一象限にある円の面積は、全体の 1/4 だから π/4。

モンテカルロ法では、第一象限の正方形の中に、ランダムに点(x,y)を打つ。たくさんのランダムな点を、疑似乱数から生成しよう。そのとき、状態を持つのではなく、乱数の無限数列を生成する。

import Random

randomSeq :: Int -> [Double]
randomSeq seed = randomRs (0,radius) (mkStdGen seed)

この乱数列を、乱数の点に変換するために、組を作る関数を定義する。

pair :: [a] -> [(a, a)]
pair (x:y:zs) = (x,y):pair zs
pair _ = error "pair"

この二つの関数から、ランダムな点の無限数列を生成する関数を作成できる。

type Point = (Double,Double)

randomPoints :: Int -> [Point]
randomPoints = pair . randomSeq

適当な初期値を与えれば、以下のように動く。

> take 5 $ randomPoints 1234
[(9.764003220424977e-3,0.26272692618955085),(0.26353785948910236,4.618864845162923e-2),(0.16897377573628297,0.8421564352098285),(0.5194425471172301,0.6740836520153293),(0.18459976992676963,0.6367203037107177)]

さてさて、ランダムに打った点が、円の中に入る確率は、面積を考えると、π/4である。点が円に入るか否かは、原点からの距離を計算すればよい。

monteTest :: Point -> Bool
monteTest (x,y) = x*x + y*y <= radius

円周率を近似的に求める関数は、ランダムな点の数と、円に収まった数の比を取ればよいから、以下のように実装できる。

calcPi :: Int -> [Point] -> Double
calcPi count ps = 4 * fromIntegral (pass ps) / fromIntegral count
  where
    pass = length . filter id . map monteTest . take count

実際に動かしてみる。

calcPi 10000 $ randomPoints 1234
→ 3.1392

というわけで、状態を持たずに疑似乱数を扱うことができた。