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

Pythonによるモンテカルロ法入門（2014/6/20）の受理・棄却法（2014/7/12）4回目。今回は練習問題2.18を解いてみた。下のような〜分布という名前がついていない完全に任意な目標分布 $f(x)$ にしたがう乱数を受理・棄却法で生成してみよう。

$\displaystyle f(x) \propto \exp (-\frac{x^2}{2}) \{ \sin (6x)^2 + 3 \cos (x)^2 \sin (4x)^2 + 1 \}$

提案分布 $g(x)$ には標準正規分布を使う。

$\displaystyle g(x) = \frac{1}{\sqrt{2 \pi}} \exp (-\frac{x^2}{2})$

目標分布の式が示すように左辺は右辺に「比例する」だけで、積分が1になるように正規化されていない。乱数を生成するだけなら正規化は必須でないが、今回は最初から正規化しておいた。正規化するには右辺を全区間で積分した値で $f(x)$ を割ればよい。Pythonで定積分を計算するには、scipy.integrate.quadが使える。

さっそく実装してみよう。

# 任意の確率分布にしたがう乱数を受理・棄却法で生成する

import numpy as np
import matplotlib.pyplot as plt
import scipy.optimize
import scipy.integrate
from scipy.stats import norm

np.random.seed()

# 正規化前の目標分布f
f = lambda x: np.exp(-x**2 / 2) * (np.sin(6*x)**2 + 3 * np.cos(x)**2 * np.sin(4*x)**2 + 1)

# 正規化定数
I = scipy.integrate.quad(f, -Inf, Inf)[0]

# 正規化済み（積分が1）の目標分布
f = lambda x: np.exp(-x**2 / 2) * (np.sin(6*x)**2 + 3 * np.cos(x)**2 * np.sin(4*x)**2 + 1) / I

# 提案分布g
gv = norm(loc=0.0, scale=1.0)
g = gv.pdf

# 分布の上限を指定する定数Mを設定
xopt = scipy.optimize.fmin(lambda x: - f(x) / g(x), 0.0, disp=False)[0]
# そこでの値をMとする
M = f(xopt) / g(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(-10, 10, 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.xlim((-5, 5))
plt.show()

実行すると

M = 1.85606973547
サンプル数: 100000 => 53818
実際の受理率  : 0.538180
理論的な受理率: 0.538773

f:id:aidiary:20140717213835j:plain

赤色が目標分布 $f(x)$ 。緑色が提案分布 $g(x)$ の正規分布。青色が受理されたサンプルのヒストグラム。ヒストグラムは赤色の理論値と一致していることがわかる。

少し気になってたのだけれど、どんな分布でも理論的な受理率が $1/M$ で計算できてしまうのが不思議に感じてきた。もう少し理論的な背景も調べた方がよさそうだな。

今回で受理・棄却法はおしまい。次回から第3章のモンテカルロ積分にいく予定。乱数生成で難しい積分が効率的に計算できてしまうのだ！