読者です 読者をやめる 読者になる 読者になる

ピュアオーディオ的な解像度で語る乱数のコク

周回遅れも甚だしいですが、「乱数のコク」について思う所があったので計算をしてみました。

オリジナル

一様乱数と比べて、5つの一様乱数の平均の方がコクがあるとのことです。 シミュレーションで得たヒストグラムは以下のようになりました。

def rand_koku(n, seed):
    np.random.seed(seed)
    u1 = np.random.rand(n)
    u2 = np.random.rand(n)
    u3 = np.random.rand(n)
    u4 = np.random.rand(n)
    u5 = np.random.rand(n)
    return (u1, (u1+u2+u3+u4+u5)/5.0)

f:id:marugari2:20161108002149p:plain

平均を取った側は正規分布のような釣り鐘型になっています。 果たしてこれをコクと呼べるでしょうか。2次モーメントも異なっていますし、分布の形も別物です。 オーディオで言ったらSONYBOSEくらい差があります。

乱数のコクが感じられる例

今回は2種類の正規乱数を用意しました。

1つはBox-Muller法による正規乱数です。 これはnumpyの正規乱数生成にも用いられているベーシックなものです。

def normal_box_muller(n, seed):
    np.random.seed(seed)
    nm = int(n / 2)
    u  = np.random.rand(n)
    u1 = u[:nm]
    u2 = u[nm:n]
    r  = -2.0 * np.log(u1)
    v  = 2.0 * np.pi * u2
    z1 = np.sqrt(r) * np.cos(v)
    z2 = np.sqrt(r) * np.sin(v)
    z  = np.concatenate([z1, z2, -z1, -z2])
    return z

もう1つはMoroの逆関数法による正規乱数です。 正規分布の累積分布関数の逆関数を利用して、一様乱数から正規乱数を得ることができます。

def normal_moro(n, seed):
    np.random.seed(seed)
    a0  = 2.50662823884
    a1  = -18.61500062529
    a2  = 41.39119773534
    a3  = -25.44106049637
    b0  = -8.47351093090
    b1  = 23.08336743743
    b2  = -21.06224101826
    b3  = 3.13082909833
    c0  = 0.3374754822726147
    c1  = 0.9761690190917186
    c2  = 0.1607979714918209
    c3  = 0.0276438810333863
    c4  = 0.0038405729373609
    c5  = 0.0003951896511919
    c6  = 0.0000321767881768
    c7  = 0.0000002888167364
    c8  = 0.0000003960315187
    u   = np.random.rand(n)
    y   = u - 0.5
    f1  = (np.abs(y) < 0.42)
    y1  = y[f1]
    r   = y1 ** 2
    num = ((a3*r+a2)*r+a1)*r+a0
    den = (((b3*r+b2)*r+b1)*r+b0)*r+1.0
    x1  = y1 * num / den
    f2  = np.logical_not(f1)
    u2  = u[f2]
    y2  = y[f2]
    r   = np.where(y2 > 0.0, 1.0 - u2, u2)
    r   = np.log(-np.log(r))
    x2  = c0+r*(c1+r*(c2+r*(c3+r*(c4+r*(c5+r*(c6+r*(c7+r*c8)))))))
    x2  = np.where(y2 < 0.0, -x2, x2)
    z1  = np.concatenate([x1, x2])
    z   = np.concatenate([z1, -z1])
    return z

アルゴリズムはどちらも Paul Glasserman, Monte Carlo Methods in Financial Engineering を参考にしています。

乱数を200,000,000個発生させて得たヒストグラムは以下のようになりました。

  • 両側0.5シグマ f:id:marugari2:20161108003322p:plain 両側0.5シグマの範囲では差を見て取ることができません。

  • 片側4シグマ f:id:marugari2:20161108003533p:plain 同じ一様乱数系列を用いていますが、正規分布への変換法によって片側4シグマの裾の分布が異なっています。

英語版WikipediaではBox-Muller法の弱点について記載があります。 裾の分布に拘る場合はMoroの逆関数のようなものを使うことが推奨されているようです。

何気なくrandnで正規乱数を得ることも多いですが、そのアルゴリズムの違いによって得られる結果が異なることを見てきました。 こういったところにこそ「乱数のコク」があるのではないでしょうか。

最後に

正規分布は非常に基礎的な確率分布ですが、それに従う乱数を得るためのアルゴリズムは(私にとっては)そんなに簡単ではありません。 正規乱数でさえそうなのです。

ベイズ統計では、複雑な事前分布と複雑な尤度関数の積からなる事後分布に従うサンプルを生成したいケースがあるのですが、そのような複雑な分布関数にはBox-Mullerのような解析的なアプローチもMoroのような精度の良い近似関数も存在しないことが常です。 そのようなマトモでない分布に従う乱数を生成するためにMarcov Chain Monte Carloが必要となっている訳です。

こういった話はRobert and Casella, Monte Carlo Statistical Methodsに感動的にわかりやすく書かれているので、私は強くオススメをしています。