乱数による正規分布の生成

(generating normally distributed random numbers)

使い方：

start/stop ボタンを押すだけです。clear ボタンで再び初期化します。
標準偏差を 50 と設定しています。標準偏差の１毎に度数を取り、横方向のピクセルとして表示しています。縦方向の１ピクセルは 10000 カウントです。
標準偏差の１倍（つまり 50）毎にも度数を取ってテキストボックスに表示しており、全体で１になるはずです。グラフの縦罫線の間隔は標準偏差毎です。

解説：

シミュレーションで乱数を利用する場合がある。モンテカルロ法などと言われるが、乱数の最も基本的なものは、有限区間、普通はゼロから１までの間の実数の均一な生成である。
乱数で次に使いたいのが、正規分布である。しかしこれはなかなか難しく、すぐに方法は見つからなかった。Wikipedia には文献が書かれているが、取り寄せて期待はずれでは残念だし。
そこで、色々考える中でできてしまったので、ご紹介する。
そうしたモンテカルロを何に使うかって？　それは続く記事を乞うご期待です。
まず正規分布から飾りを除いた基本部分は次の関数である。
　　　

これは確率密度関数ではない。確率密度関数とは、全領域での積分が１となるものと定義される。また、標準偏差というのを好みに合わせる必要があるが、とりあえず、これをいじってみる。

まず、積分値を知りたいが、意外にもこれは積分の２乗をいきなり求めることによって得られる。
　　　

従って、
　　　

次に、標準偏差を得るための積分をしたいが、折角なので一般論として次の積分を求めておきたい。
　　　

先に求めたものはこの最初の場合である。
　　　

次のインスタンスはさらに簡単である。
　　　

n=3 からのものは、それぞれ n-2 の結果をもとに、部分積分によって再帰的に求めることができる。
　　　

さて、さきほど

を求めるにあたって、２次元の極座標からその平方を求めたわけだが、一般に

の n 乗を n 次元の極座標から計算する手順を考える。今度は

を書き換えることはせず、極座標に変換するだけである。超球の表面積を
　　　

と置くと、
　　　

この最後の積分を

の定義と見比べると、その半分であることがわかる。そこで、
　　　

整理すると次の表のようになる。

は n/2 を引数とするガンマ関数に等しい。

表：次元ｎと積分値

さて、表の n=1 のバリエーションとして、分散が σ^2 の確率密度関数（積分値＝１）を求めるなら、
　　　

これを解くと、
　　　

確率密度関数は、
　　　

となる。
なお、熱力学との対応は、

すなわち、速度成分と置くと、
　　　

となる。

これを乱数で生成しようとすると、一般には、確率密度関数を積分した累積分布関数というのを求める。これは、変位から、累積頻度を求めるもので、この累積分布関数に対する逆関数（頻度から変位を求める）に、一様分布の乱数を放り込むと、各回の値が求められるという理屈になる。正規分布の場合のこの「逆関数」は Probit と言うらしいが、これはとても難しい手続きのように思われる。
しかし、先ほど積分を、一度２次元化して、極座標で求めたような方法で解決できるのではないかと考えた。言わば、一次元の Probit を求めなくても、極座標化した二次元の逆関数から二次元の正規分布は作ることができ、それを一軸に射影すれば、一次元の乱数になるということである。

乱数生成部分のプログラムは次のような簡単なものである（最初は簡単だったが、一様分布部分の精度が少し悪く、改善すると共に複雑になった）。

また、呼び出す側は次のような感じ。

const distr=new rannormal(50);
...
var x=distr.step();