Box-Muller法

Pythonによるモンテカルロ法入門（2014/6/20）のつづき。今回はBox-Muller法を使って標準正規分布に従う乱数を生成してみました。あとおまけで中心極限定理に基づく正規乱数の生成も試してみます。

Box-Muller法は、2つの一様乱数から標準正規分布 N(0, 1) に従う2つの乱数を生成するアルゴリズム。 $U_1$ と $U_2$ が独立同分布の $U_{[0,1]}$ とすると、次の2つの変数

$X_1 = \sqrt{-2 \log (U_1)} \cos (2 \pi U_2)$
$X_2 = \sqrt{-2 \log (U_1)} \sin (2 \pi U_2)$

は独立同分布の $N_{[0,1]}$ になるらしい。さっそく検証してみよう。

# Box-Mullerアルゴリズムを用いて標準正規分布N(0,1)に従う乱数を生成
import numpy as np
from scipy.stats import uniform, norm
import matplotlib.pyplot as plt

# 独立した一様分布からそれぞれ一様乱数を生成
np.random.seed()
N = 100000
rv1 = uniform(loc=0.0, scale=1.0)
rv2 = uniform(loc=0.0, scale=1.0)
U1 = rv1.rvs(N)
U2 = rv2.rvs(N)

# Box-Mullerアルゴリズムで正規分布に従う乱数に変換
# 2つの一様分布から2つの標準正規分布が得られる
X1 = np.sqrt(-2 * np.log(U1)) * np.cos(2 * np.pi * U2)
X2 = np.sqrt(-2 * np.log(U1)) * np.sin(2 * np.pi * U2)

# 変換した乱数の分布と標準正規分布の真のpdfを描画
rv = norm(loc=0, scale=1)

# X1の分布
plt.subplot(2, 1, 1)
nbins = 50
plt.hist(X1, nbins, normed=True)
x = np.linspace(rv.ppf(0.01), rv.ppf(0.99), 1000)
y = rv.pdf(x)
plt.plot(x, y, 'r-', lw=2)
plt.xlim((rv.ppf(0.01), rv.ppf(0.99)))

# X2の分布
plt.subplot(2, 1, 2)
plt.hist(X2, nbins, normed=True)
x = np.linspace(rv.ppf(0.01), rv.ppf(0.99), 1000)
y = rv.pdf(x)
plt.plot(x, y, 'r-', lw=2)
plt.xlim((rv.ppf(0.01), rv.ppf(0.99)))

plt.show()

実行すると

f:id:aidiary:20140706194346p:plain

確かにX1とX2の乱数は標準正規分布と一致していることが確認できました。

中心極限定理に基づく近似アルゴリズム

もう一つ紹介されていた中心極限定理（CLT: Central Limit Theorem）に基づく方法も試してみました。中心極限定理に基づく方法はBox-Muller法に比べて精度は落ちるけど高速とのこと。

練習問題2.3と合わせて[-0.5, 0.5]の一様乱数12個を足し合わせることで1つの正規乱数を生成してみます。実際には正規乱数を1つ生成するのではなく、100000個まとめて生成してヒストグラムを描いています。

# 中心極限定理（CLT）による正規乱数生成器
import numpy as np
from scipy.stats import uniform, norm
import matplotlib.pyplot as plt

# 独立した一様分布からそれぞれ一様乱数を生成
np.random.seed()
N = 100000

# 12個の一様乱数を使用
K = 12

# (-0.5, 0.5)の一様乱数を生成
rv = uniform(loc=-0.5, scale=1.0)
U = rv.rvs(K * N).reshape((K, -1))

# K個の一様分布の和から正規分布を生成
Z = np.sum(U, axis=0)

# 変換した乱数の分布と標準正規分布の真のpdfを描画
plt.figure(1)
nbins = 50
plt.hist(Z, nbins, normed=True)

rv = norm(loc=0, scale=1)
x = np.linspace(rv.ppf(0.01), rv.ppf(0.99), 1000)
y = rv.pdf(x)
plt.plot(x, y, 'r-', lw=2)
plt.xlim((rv.ppf(0.01), rv.ppf(0.99)))

plt.show()