人工知能に関する断創録

このブログでは人工知能のさまざまな分野について調査したことをまとめています(更新停止: 2019年12月31日)

受理・棄却法 (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章のモンテカルロ積分にいく予定。乱数生成で難しい積分が効率的に計算できてしまうのだ!