曲がった迷路その16 - 一様でない乱数の発生

前回離散的な2次元ポテンシャルの補間を考えた。その中にある分布を持った長さの壁を作る。その分布をどうするかを考える。

壁の長さの分布

今回のアルゴリズムでは乱数を使う場所は2カ所ある。ひとつはスペックルのもとになる反射率を作るところ、もうひとつはひとつの壁の長さを決めるところ。反射率は一様な乱数がふさわしいけど、壁の長さは一様な分布だとランダム性が大きすぎてやはり無機的な感じになる。こういうときは正規分布のようにある値のまわりに広がった分布が望ましい。

そういう一様でない乱数を生成する方法をちょっと考えてみる。

一様でない乱数の生成

疑似乱数の生成とその品質

普通、疑似乱数は任意の整数を一様に発生させるものが多い。例えばunixの古典的なrand()関数は0からRAND_MAXまでの自然数をほぼ一様に発生させる。rand()関数は整数に対する非線形な操作を繰り返し作用させることで整数列を発生させて乱数としている。「非線形な操作」は決定論的（たいてい割った余りを使う）なもので、最初の整数（seedと呼ぶ）を決めると同じ整数列が得られることになる。従って厳密には「疑似乱数」と呼ばれる。

疑似乱数の品質は

• 一様であること
• 前後が無相関であること
• 生成効率

などで評価される。unixの場合、もとのrand()関数の品質が問題になり、いくつかの改良版が存在している。特に最近は暗号化の絡みで高い品質の乱数発生器が要求されることがあるが、そういった応用にはunixのrand()系の疑似乱数は向かない。

また、乱数を整数のまま使う用途では偶数奇数の分布やある桁の0と1の分布などに偏りが問題になることがある。そのような用途でもrand()系の関数は使えない。

今回のような応用では、一様な実数が必要となる。

rand()関数のような自然数の発生器を使って区間[0,1)の浮動小数点を発生させるには

double frand()
{
    return (double)rand() / RAND_MAX;
}

みたいにしてやるのが簡単。これだと、[0,1)がだいたい一様（とびとびだけど）に現れることになる。

unixにはrand()を改良したrandom()や区間[0,1)のdoubleを返すdrand48()という関数もある。が、世間的な評価ではどの発生器も五十歩百歩のようで、man pageにはどれも使わない方がよいなどと書かれていることがある。特にrand()の古い実装では偶数と奇数がかわりばんこに現れるようなものもあるらしい。

今回のような場合には品質はあまり問題にならないが、気にしておいた方がよい。

一様でない乱数

一様ではなくて例えば発生確率が正規分布（ガウシアン）になるようにするにはどうすればいいか。

gslなどのライブラリやMathematica、MATLABなどにははじめから備わっているけど、この問題に出会うたびに悩むのでここに書いておく。

といっても、順を追えばそれほど難しくない。
値域が区間I=[a,b]の乱数で、ある値x∈Iとそのすぐ近くx+δの間の数が発生する確率Pが

となるようなものR_fをI₀=[0,1]区間の一様な乱数の発生器R_Uから作ることを考える。

ここで

である。このf(x)を確率密度関数とよぶ。

まず、

というような関数g(z)を考える。これはz∈[0,1]から[a,b]への関数で、一様乱数の実数を入力すると、その分布が式-42となるようなもので、

のようなものである。
図-23にイメージを示す。

横軸が入力となる一様乱数の区間でg(z)の関数によって、縦軸の分布f(x)の密度になる。

この図からわかるように幅がf(x)の逆数に変換されればいい。ということは

であればいいということになる。

式-45はわかりにくいのでg(z)の逆関数

を考える。

図から式-45を考えたときと同じようにする（縦横を取り替える）と

であればよくて、これはすぐ解けて

となる。ここでc₁とc₂は比例定数と積分定数で、従って

であればいい。 つまり、分布関数を積分して逆関数を求めればいい、ということになる。そしてc₁とc₂は

となるように決めれば良い。