プログラムな日常>ウェブプログラム作品集>
乱数による正規分布の生成(2024)
(generating normally distributed random numbers)
 
sigma_result/sigma_set percent/count

2年前に作ったページの改訂版です。見た目の変化は少ないが、「渋い」改訂してます。 主な改訂点は、


使い方:
解説:

シミュレーションで乱数を利用する場合がある。モンテカルロ法などと言われるが、乱数の最も基本的なものは、有限区間、普通はゼロから1までの間の実数の均一な生成である。

乱数で次に使いたいのが、正規分布である。色々考える中で、解法と証明ができたので、ご紹介する。

まず正規分布から飾りを除いた基本部分は次の関数である。
   
これは確率密度関数ではない。確率密度関数とは、全領域での積分が1となるものと定義される。また、標準偏差というのを好みに合わせる必要があるが、とりあえず、これをいじってみる。

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

このように、全空間についての積分が2次元とすることで簡単に求められたが、2次元でなら半径に沿って一定範囲の積分も求めやすいのが特筆すべき点である。

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

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

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

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

さて、さきほどB1を求めるにあたって、2次元の極座標からその平方を求めたわけだが、一般に B1の n 乗を n 次元の極座標から計算する手順を考える。 超球の表面積を
   
と置き、極座標に変換する。
   

この最後の辺の積分部分をBnの定義と見比べると、正負で対称な関数につき積分範囲が半分であることがわかる。そこで、
   

整理すると次の表のようになる。Bn は n/2 を引数とするガンマ関数に等しい。


表中で、変数Anは超球の定義
   
に使われる係数、Bnは次の関係式で与えられる。
   

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

これを解くと、
   

確率密度関数は、
   
となる。

なお、熱力学との対応は、x=v すなわち、速度成分と置くと、
   
となる。

これを乱数で生成しようとすると、一般には、確率密度関数を積分した累積分布関数というのを求める。 これは、変位から、累積頻度を求めるもので、この累積分布関数に対する逆関数(頻度から変位を求める)に、一様分布の乱数を放り込むと、各回の値が求められるという理屈になる。 正規分布の場合のこの「逆関数」は Probit と言うらしいが、これはとても難しい手続きのように思われる。

しかし、先ほど全空間の積分を、一度二次元化して、極座標で求めたように、一次元の Probit を求めなくても、極座標化した二次元の逆関数から二次元の正規分布を作り、それを一軸に射影することで、一次元の乱数にできる。

乱数生成部分のプログラム("../_common/ran.js")は次の通り。

描画のプログラム("001.js")は次の通り。

その中で、乱数(正規分布)に対するコールは次の通り。
const distr=new rannormal();
...
var x=distr.step()*sigma;

C#版の乱数生成プログラム("Ran.cs")は次の通り。