JavaScriptでBox-Muller法による正規分布からのサンプリング

1. 背景

標準正規分布に従う確率変数を生成する方法としてBox-Muller法（ボックス＝ミュラー法）が知られている。
Box-Muller法を用いると、一様分布に従う確率変数を変換することで正規分布に従う擬似乱数を発生させることができる。一様分布に従う乱数は殆どのプログラミング言語で提供されている（JavaScriptならMath.random()）ため、Box-Muller法と組み合わせれば正規分布からのサンプリングが可能になる。

2. 数式

$U_1, U_2$が区間$(0, 1)$における一様分布に従う互いに独立な確率変数のとき、以下の式によって定義される$Z_0, Z_1$は、標準正規分布$N(0, 1)$に従う互いに独立な確率変数となる。

\left\{
    \begin{array}{l}
     Z_0 = R cos(\Theta) = \sqrt{-2lnU_1} cos (2 \pi U_2) \\
     Z_1 = R sin(\Theta) = \sqrt{-2lnU_1} sin (2 \pi U_2)
    \end{array}
  \right.

3. 実装

標準で利用できるMath.random()が一様分布に従う乱数であることを利用する。

// Generate n samples from standard normal distribution
function boxMuller(n) {
  const samples = [];
  Array(Math.ceil(n / 2)).fill().forEach(_ => {
    const R = Math.sqrt(-2 * Math.log(Math.random()));
    const theta = 2 * Math.PI * Math.random();
    samples.push(R * Math.cos(theta)); // z1
    samples.push(R * Math.sin(theta)); // z2
  });
  // if n is odd, drop the last element
  return samples.slice(0, n);
}

例えばboxMuller(10)を呼ぶと、標準正規分布からのサンプルが10個収められたArrayが返ってくる。

4. デモ

実装はObservableというサイトのノートブック: Box-Muller transformに公開しているため、ブラウザ上で弄ることが可能。

上記関数によって生成された10000個のサンプルについて、

• ヒストグラム
• 正規Q-Qプロット

を描画したものを以下に示す。

ヒストグラム

正規Q-Qプロット
ほぼ直線上にあることが分かる。