受理・棄却法 (2) - 人工知能に関する断創録

Pythonによるモンテカルロ法入門（2014/6/20）の受理・棄却法（2014/7/12）の2回目。前回は、提案分布に一様分布を使ったためサンプリングの効率が悪かった。今回は、より効率のよい提案分布gとしてベータ分布を使ってみよう。目標分布fは前回と同じベータ分布。下の図を見ればわかるように提案分布では目標分布を覆えていない。すっぽり覆えるように提案分布をM倍する必要がある。

f:id:aidiary:20140713213841p:plain

Mを求めるには前回と同様にscipy.optimizeを使った。

# f(x)/g(x)を最大化するxoptを求める
xopt = scipy.optimize.fmin(lambda x: - f(x) / g(x), 0.0, disp=False)
# そこでの値をMとする
M = f(xopt) / g(xopt)
print "xopt =", xopt
print "M =", M

前回は g(x) が一様分布でつねに1だったので f(x) の最大化だけでよかった。今回は、g(x) が定数ではないので f(x)/g(x) の比を最大化する x を求める必要がある。結局、x = 0.7 で f(x) / g(x) の比が最大になっており、M = 1.67 とすればすっぽり覆えることがわかった。受理・棄却法で目標分布にしたがう乱数を生成してみる。

# ベータ乱数を受理・棄却法で生成
# 目標分布（ここではベータ分布）のpdfは既知とする
# 提案分布として目標のベータ分布を覆うベータ分布を使用

import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize
from scipy.stats import uniform, beta

np.random.seed()

# 目標分布f
f = beta(a=2.7, b=6.3).pdf

# 提案分布g
gv = beta(a=2.0, b=6.0)
g = gv.pdf

# 分布の上限を指定する定数Mを設定
# f(x)/g(x) <= Mを満たす必要がある

# f(x)/g(x)を最大化するxoptを求める
xopt = scipy.optimize.fmin(lambda x: - f(x) / g(x), 0.0, disp=False)
# そこでの値をMとする
M = f(xopt) / g(xopt)
print "xopt =", xopt
print "M =", M

# 受理・棄却法
Nsim = 100000

# 提案分布gからの乱数Yを生成
Y = gv.rvs(size=Nsim)

# 一様乱数UをNsim個生成
U = uniform.rvs(size=Nsim)

# Yから受理の条件を満たすサンプルXを残して残りを棄却
X = Y[U <= f(Y) / (M * g(Y))]
print u"サンプル数: %d => %d" % (len(Y), len(X))
print u"実際の受理率  : %f" % (len(X) / float(len(Y)))
print u"理論的な受理率: %f" % (1.0 / M)

# 目標分布を描画
x = np.linspace(0.0, 1.0, 1000)
y = f(x)
plt.plot(x, y, 'r-', lw=2)

# 提案分布（ベータ分布）を描画
y = M * g(x)
plt.plot(x, y, 'g-', lw=2)

# 受理した乱数の分布を描画
plt.hist(X, bins=50, normed=True)

plt.show()

実行すると

xopt = [ 0.7]
M = [ 1.67180777]
サンプル数: 100000 => 59817
実際の受理率  : 0.598170
理論的な受理率: 0.598155

f:id:aidiary:20140713213843p:plain

ベータ分布にしたがう乱数が生成できていることがわかる。10万サンプル生成して59817サンプルが受理されているので受理率は59%となる。理論値の1/Mとも一致している。前回、提案分布に一様分布を使ったときの
受理率は37%だったので少し改善している。

実際に生成したサンプルを図示してみよう。

# 受理されたサンプルと棄却されたサンプルを描く
# 点の数が多すぎるのでNsimを小さくした

# 目標分布f
f = beta(a=2.7, b=6.3).pdf

# 提案分布g
gv = beta(a=2.0, b=6.0)
g = gv.pdf

Nsim = 2000

# 候補密度からの乱数Yを生成
Y = uniform.rvs(size=Nsim)

# 一様乱数UをNsim個生成
U = uniform.rvs(size=Nsim)

# 受理されたサンプルと、棄却されたサンプルのインデックスを計算
acceptedIdx = U <= f(Y) / (M * g(Y))
rejectedIdx = U > f(Y) / (M * g(Y))

# 目標分布を描画
x = np.linspace(0.0, 1.0, 1000)
y = f(x)
plt.plot(x, y, 'r-', lw=2)

# 提案分布を描画
y = M * g(x)
plt.plot(x, y, 'g-', lw=2)

# 受理されたサンプルを描画
plt.scatter(Y[acceptedIdx], U[acceptedIdx] * M * g(Y[acceptedIdx]), color="red")
plt.scatter(Y[rejectedIdx], U[rejectedIdx] * M * g(Y[rejectedIdx]), color="blue")

plt.xlim((0.0, 1.0))
plt.ylim((0.0, 5.0))

plt.show()