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
というわけで、状態を持たずに疑似乱数を扱うことができた。