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

前回、一様乱数から指定された出現確率に従った乱数を発生させるための数学を考えた。今回、それにもとづいて普段よくやる例ふたつをみる。

単位円内の例

よくあるのが、2次元で、単位円内に一様に乱数を発生させたい場合。

簡単な方法は

として

を満たす組みだけ取り出す、というものである。ここでR_Uは昨日出てきた[0,1]区間の一様乱数発生器。簡単だけど、これでは5回に1回は捨てることになって効率が悪い。品質のいい乱数は発生コストも高いのであまり捨てたくない。
そこでうっかり

としてしまうが、当然これはだめ。中心付近の出現確率が高くなってしまう。

これは、ちゃんとヤコビアンに従って

で、R_Lの発生確率f_L(r)が

としなければいけない。
これを一様乱数から作る変換関数g_L(z)は式-49から、

となる。

ガウシアンの例

ガウシアン、いわゆる正規分布に乱数を発生させることを考える。 この場合、確率分布関数f_G(x)はaを正の定数として

である。この変換関数g_G(z)は簡単には求まらない。無理矢理書くと

となる。ここでerf^-1(x)は誤差関数の逆関数で、誤差関数erf(x)とは

である。そのまんま。なんにもやってないのと同じだけど、みんなよく使うので数値評価はされている。MathematicaやMATLABは標準関数として数値評価できる。Mathematicaはどうやら誤差関数を数値的に解き直しているらしい。

実用的なアルゴリズムとしてボックス・ミューラー法という方法がある。これはいろいろなところに解説があるので省略する。gslもガウシアンの乱数発生器はボックスミューラーを使っている。

他に、ベキ展開してしまう、と言う方法もある。
生の逆誤差関数erf^-1(x)をMathematicaでx=0で10次までのベキ展開すると

となる。x=0の近くでは素直な形をしているので残差は小さい。

x=1で展開するとMathematicaは面白いものを返してくる。

なんでこんなものが出てくるのかわからないが、この二つの式をx=0.8あたりでつなげば数値的には残差は10^-3程度に収まって、乱数発生器とするには十分な精度が得られる。もっとこまかく区間に分ければ次数も低くて高い精度が得られる。もちろんボックス・ミューラーの方が効率はいいけど、こういうのもたまには面白い。

特殊な分布（確率密度関数）の乱数を発生させたい場合などで、積分や逆関数が解析的に求まらない場合、こういう力づくの手もある。