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

Pythonによるモンテカルロ法入門（2014/6/20）の受理・棄却法（2014/7/12）の3回目。今回は、練習問題2.8を解いてみた。目標分布を標準正規分布、提案分布を二重指数分布（ラプラス分布）とし、受理・棄却法（棄却サンプリング）を用いて標準正規分布にしたがう乱数を生成せよというのが課題。二重指数分布（double exponential distribution）という分布があることを始めて知った。ラプラス分布（Laplace distribution）の別名とのこと。
http://en.wikipedia.org/wiki/Laplace_distribution

目標分布f（標準正規分布）と提案分布g（ラプラス分布）は下の式で与えられる。

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

$\displaystyle g(x|\alpha) = \frac{\alpha}{2} \exp (-\alpha |x|)$

今まで書いてきたように受理・棄却法を使うには目標分布をすっぽり覆うように提案分布をM倍する必要がある。このMを求めるには

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

を最大化するxを求めればよい。前回はscipy.optimizeを使って数値計算で求めたけど今回は解析的に求めてみた。まず、xで微分して0と置くと、 $x=\alpha$ で $f(x)/g(x)$ が最大値を取ることがわかる。 $x=\alpha$ を代入した値が最大値として上限を与えるので

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

が成り立つ。次にMはなるべく小さい方が棄却率が低くなり効率がよいため上限の値を最小化するパラメータ $\alpha$ を求める。 $\alpha$ で微分して0とおくと $\alpha = 1$ で最小値を取ることがわかる。 $\alpha = -1$ でも極値を取るが、提案分布のgが負の値になるため適さない。

本当にこれで合ってるのか疑問なときはグラフを書いてみるとよい。matplotlibでも描けるけどちょっとした検算はWolfram Alphaがすごく便利。検索窓に「sqrt(2/pi) a^-1 exp(a^2/2)」を入れると下のようなグラフを返してくれる。確かにa=1で最小値を取ることがすぐわかる。

これでMが求まったためさっそく受理・棄却法でサンプリングしてみよう。

# 正規乱数を受理・棄却法で生成
# 提案分布としてラプラス分布（二重指数分布）を使用

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

np.random.seed()

# 目標分布f
f = norm(loc=0.0, scale=1.0).pdf

# 提案分布g
# scale = 1/alphaになる
# 練習問題2.8からもっとも効率がよいのがalpha=1なのでscale=1
gv = laplace(loc=0.0, scale=1.0)
g = gv.pdf

# 分布の上限を指定する定数Mを設定
# f(x)/g(x) <= Mを満たす必要がある
# 練習問題2.8から下式がもっとも効率のよいMとなる
alpha = 1.0
M = np.sqrt(2.0 / np.pi) * (1.0/alpha) * np.exp(alpha**2/2.0)
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.show()

実行すると

M = 1.31548924696
サンプル数: 100000 => 76170
実際の受理率  : 0.761700
理論的な受理率: 0.760173

f:id:aidiary:20140716215209j:plain

ラプラス分布（緑色の線）を提案分布として正規乱数（青色のヒストグラム）が生成されることが確認できた。正規乱数のヒストグラムは理想的な正規分布（赤色の線）と一致している。