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

人工知能に関する断創録

人工知能、認知科学、心理学、ロボティクス、生物学などに興味を持っています。このブログでは人工知能のさまざまな分野について調査したことをまとめています。最近は、機械学習・Deep Learningに関する記事が多いです。



Box-Muller法

機械学習 PRML

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

Box-Muller法

Box-Muller法は、2つの一様乱数から標準正規分布 N(0, 1) に従う2つの乱数を生成するアルゴリズム。U_1U_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()

f:id:aidiary:20140706211747p:plain

ちゃんと正規分布に従う乱数になっていました。

補足

一様乱数を12個足し合わせると何で標準正規分布に従う乱数になるのかすぐに理解できなかったので少し補足。


中心極限定理
によると平均が\mu、分散が\sigma^2で独立同分布(i.i.d)のn個の確率変数 X_1, X_2, \cdots, X_n の平均

\displaystyle S_n = \frac{X_1 + X_2 + \cdots + X_n}{n}

は、nが大きい場合、正規分布 N(\mu, \sigma^2/n) に近づく。今回は平均\mu=0、分散\sigma^2=1/12になる一様分布を使っていて確率変数の数はn=12なので、中心極限定理より

\displaystyle S_{12} = \frac{1}{12} \sum_{i=1}^{12} X_i

は、N(0, 1/12^2)に従う。つまり、

\displaystyle 12 S_{12} = \sum_{i=1}^{12} X_i (12個の一様分布の単なる和)

は、N(0, 1)の標準正規分布に従う(証明終了)

最後に使ったのは分散の定理 V(aX) = a^2 V(X) です。12個使ったのは一様分布の分散が1/12だったからなのか・・・納得。