Qiita Teams that are logged in
You are not logged in to any team

Log in to Qiita Team
Community
OrganizationEventAdvent CalendarQiitadon (β)
Service
Qiita JobsQiita ZineQiita Blog
1
Help us understand the problem. What are the problem?

More than 1 year has passed since last update.

@sw1227

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プロット

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

ヒストグラム
スクリーンショット 2019-08-25 15.11.37.png

正規Q-Qプロット
ほぼ直線上にあることが分かる。
スクリーンショット 2019-08-25 15.11.44.png

Why not register and get more from Qiita?
  1. We will deliver articles that match you
    By following users and tags, you can catch up information on technical fields that you are interested in as a whole
  2. you can read useful information later efficiently
    By "stocking" the articles you like, you can search right away
1
Help us understand the problem. What are the problem?