F#でクロージャを使ったBox-Muller法を書いた

Box-Muller法で正規分布$\mathrm{N}(0, 1)$に属する乱数を生成するF#コードを書きました。
標準の.NETには$\mathrm{U}[0, 1)$の疑似一様乱数を返す関数はありますが、正規分布を返す関数はありません。
数値計算ライブラリには当然入っているのですが、正規分布のためだけに導入するのは大げさなので自分で書きました。

F#の勉強も兼ねて、クロージャを使って書いてみました。

Box-Muller法

以下の式を用いると、独立な一様分布乱数（以下、一様乱数）$a, b \in \mathrm{U}(0, 1)$から独立な正規分布乱数（以下、正規乱数）$x, y \in \mathrm{N}(0, 1)$を生成することができます。

x = \sqrt{-2 \log a} \cos(2 \pi b)\\
y = \sqrt{-2 \log a} \sin(2 \pi b)

クロージャを使う目的

ここでは2つの目的でクロージャを使います。
1つは余分に生成される正規乱数を保持するため、もう1つは乱数ジェネレータを隠蔽するためです。

余分に正規乱数が生成されるのは、Box-Muller法が一度に乱数を2つ生成するためです。
Box-Muller法の式を単純に実装した関数は呼び出しの度に2つの正規乱数$x, y$を生成します。
しかし通常の使い方では、正規乱数を1つだけ使うために関数を呼び出すことが多いと思います。
仮に乱数$x$を使った場合、使われなかった$y$は捨てられることになります。

なんだかもったいないと思いませんか？
え、思わない……？

もったいない以外の理由として、片側の乱数を捨てることができる実装は疑似乱数の再現性を崩す場合があります。
例えば次の疑似コードで関数simA, simBはどちらも同じ初期状態の正規乱数生成器から乱数を6個取得する処理を書いていますが、返ってくる数列は異なります。

// newNormalRandomは受け取ったシード値を元に、「呼び出しごとに2つの正規乱数を返す関数」を生成する関数とします。
let newNormalRandom : seed:int -> (unit -> float * float) = (* 実装していないので動きません。 *)

let simA () =
    let normal : unit -> float * float = newNormalRandom 42
    let normalSeries = Array.zeroCreate 6
    for i in 0 .. 5 do
        let (x, _) = normal ()
        normalSeries.[i] <- x
    normalSeries

let simB () =
    let normal : unit -> float * float = newNormalRandom 42
    let normalSeries = Array.zeroCreate 6
    for i in 0 .. 2 do
        let (x, y) = normal ()
        normalSeries.[2*i] <- x
        normalSeries.[2*i+1] <- x
    normalSeries

assert (simA () <> simB ())

この差は捨てられる乱数の扱いが異なるからです。
使い方に気をつけるという手もありますが、そもそも値を捨てれない実装にするほうが安全です。

使われなかった乱数を保持して次の呼び出しに使いまわす場合、乱数生成に使った乱数ジェネレータが外部から触られないことを保証しなければなりません。
使われなかった乱数を保持するということは、同時に生成された乱数$x, y$について$x$を利用してから次に$y$を使うまで乱数ジェネレータの状態が変化しないことを暗に仮定しています。
この仮定が満たされない場合、次の疑似コードのような、乱数の生成順序が異なるにも関わらず同じ乱数が生成されるコードが書けてしまいます。

// normalは乱数ジェネレータrngを使って生成した1つの正規乱数を返す関数とします。
// しかしながら実装は、奇数回の呼び出しで乱数を2つ生成して1つを返し、
// 偶数回の呼び出しでは奇数回で使わずに保持しておいた乱数を返しているとします。
let normal (rng : Random) : float = (* 実装していないので動かない。 *)

let simC () =
    let rng = new Random(42)
    let x = normal rng
    let y = normal rng
    let u = rng.NextDouble()
    (x, y, u)

let simD () =
    let rng = new Random(42)
    let x = normal rng
    let u = rng.NextDouble()
    let y = normal rng
    (x, y, u)

assert (simC () = simD ())

つまり乱数ジェネレータを外部におく実装も疑似乱数の再現性を崩す場合があります。

そこでクロージャを使い、余った乱数と乱数ジェネレータのスコープを制限します。

実装

uniform関数はRandom.NextDoubleメソッドをラップし、乱数が0の場合に再抽選を行います。
そのほかは上で説明した通りです。

open System

// 一様乱数U(0, 1)を返す。
let rec uniform (rng : Random) : float =
    match rng.NextDouble() with
    | 0.0 -> uniform rng
    | r -> r

// 正規乱数を返すクロージャを生成する。
let generateNormalGenerator (seed : int) : (unit -> float) =
    // 乱数ジェネレータ。
    let rng = new Random(seed)
    // 余分に生成された正規乱数を次回呼び出しまで保持するメモリ。
    let mutable memory : float option = None
    (fun () ->
        // メモリに正規乱数が保持されているかで分岐。
        match memory with
        | Some y ->
            // メモリからポップ。
            memory <- None
            y
        | None ->
            // Box-Muller法で正規乱数を生成。
            let (a, b) = (uniform rng, uniform rng)
            let x = Math.Sqrt(-2. * Math.Log a) * Math.Cos(2. * Math.PI * b)
            let y = Math.Sqrt(-2. * Math.Log a) * Math.Sin(2. * Math.PI * b)
            // 片側をメモリにプッシュ。
            memory <- Some y
            x)

次のように使います。

let normal = generateNormalGenerator 42
let size = 3000
let series = Array.init size (ignore >> normal)
let mean = Array.sum series / float size
let variance =
    let sqSeries = Array.map (fun s -> (s - mean) ** 2.0) series
    Array.sum sqSeries / float size
printfn "The first five elements: %A" series.[0 .. 4]
printfn "Mean: %f" mean
printfn "Variance: %f" variance

するとそれっぽい結果がでます。

The first five elements: [|0.5685275455; 0.6952639111; -2.01649972; -0.2904074307; -0.1491822656|]
Mean: -0.007693
Variance: 1.045807