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

Pythonによるモンテカルロ法入門（2014/6/20）のつづき。今回は、2.3節の受理・棄却法（acceptance-rejection method）、別名、棄却サンプリング（rejection sampling）による乱数生成を実験してみる。この方法を使うと事実上あらゆる分布（特定の名前がついていない分布を含めて）から乱数を生成できるようになる。

受理・棄却法では目標分布fと提案分布g*1という二つの分布を使う。目標分布fにしたがう乱数を生成するのが最終目的だが、目標分布fの乱数は直接生成できないため代わりに提案分布gの乱数を生成し、ある条件を満たす場合だけ生成した乱数を目標分布fの乱数とみなすというアプローチが受理・棄却法。条件を満たさない乱数は単に捨てられる。

このとき、目標分布fと提案分布gの関係にはいくつか仮定がある。

目標分布fにしたがう乱数は直接生成できない
逆に提案分布gにしたがう乱数は生成できる
目標分布fと提案分布gの確率密度関数は既知
すべての x について f(x) <= M * g(x) となる M が存在する

最後の条件は少しわかりにくいが、目標分布が提案分布をM倍した分布にすっぽり覆われていることを意味する。あとで挙げる図を見るとはっきりするかも。

アルゴリズム：受理・棄却法
以下の手順で生成した乱数 X は目標分布fにしたがう乱数となる。

提案分布gにしたがう乱数 Y を生成する

[0,1] の一様乱数Uを生成する

U*M*g(Y) <= f(Y)であれば Y を X として受理

そうでなければ Y を棄却し、1に戻る

Mはどうやって求めるんだ？何でこれでOKなんだ？とよくわからなかったので具体例をもとに検証してみたい。テキストには、

目標分布fをベータ分布
提案分布gを一様分布

とし、ベータ乱数を生成する例があったのでこれを実装してみた。本来、ベータ乱数は一般変換法（2014/6/26）で指数乱数から生成できるが今回は練習だから。

目標分布をすっぽり覆うMの求め方

f:id:aidiary:20140712203231p:plain

上の図に示すような目標のベータ分布fにしたがう乱数を生成したい。提案分布gには一様分布を使う。このとき、目標のベータ分布fをすっぽり覆うような f(x) <= M * g(x) を満たす M を求める必要がある。提案分布gが区間が [0,1] の一様分布ならばその最大値は1.0なので、単に目標分布 f(x) の最大値の値を M とすれば M*g(x) は f(x) をすっぽり覆えることになる。Pythonで関数 f(x) を最大化する x を求めるには最適化モジュール scipy.optimize が使える。scipy.optimize には、f(x) を最小化する x を求める fmin() しかないが、-f(x) のようにマイナスをつければ f(x) を最大化する x が求められる。この x を元の関数fに代入すれば f(x) の最大値が求められる。

# 分布の上限を指定する定数Mを設定# ベータ分布のpdfの上限値を指定すればベータ分布をすべて覆える
# 最大値を求めるためにベータ分布のpdfにマイナスをつけて
# 最小値問題に帰着させる
xopt = scipy.optimize.fmin(lambda x: -f(x), 0.0, disp=False)
M = f(xopt)[0]
print "M =", M

今回の例では、M=2.6697が得られる。

受理・棄却法の実装

先の受理・棄却法をPythonで実装してみる。

# ベータ乱数を受理・棄却法で生成
# 目標分布（ここではベータ分布）の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も保持
gv = uniform
g = gv.pdf

# 分布の上限を指定する定数Mを設定
# ベータ分布のpdfの上限値を指定すればベータ分布をすべて覆える
# 最大値を求めるためにベータ分布のpdfにマイナスをつけて
# 最小値問題に帰着させる
xopt = scipy.optimize.fmin(lambda x: -f(x), 0.0, disp=False)
M = f(xopt)[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(0.0, 1.0, 1000)
y = f(x)
plt.plot(x, y, 'r-', lw=2)

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

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

plt.show()

実行すると

M = 2.66974399495
サンプル数: 100000 => 37427
実際の受理率  : 0.374270
理論的な受理率: 0.374568

f:id:aidiary:20140712203245p:plain

受理されたサンプルのヒストグラムを描いてみると理論的なベータ分布と一致していることが確認できた。この例では、10万サンプル生成しているが、受理されるのは37427サンプルで受理率は37%にすぎない。理論的な受理率は1/Mで計算できるらしいが、それとも一致していることが確認できる。提案分布に一様分布を用いると大部分のサンプルは棄却されてしまい無駄が多いことがわかる。次回は受理率がより高い別の提案分布を使う場合を実装してみたい。

なぜこの方法でよいのか？

アルゴリズムを見ただけでなぜこれで目標分布に従う乱数が生成できるのかいまいち納得できなかったので手順を図示してみた。

(1) 提案分布gにしたがう乱数 Y を生成する

f:id:aidiary:20140712203248p:plain

(2) [0,1] の一様乱数Uを生成する
(3) U*M*g(Y) <= f(Y)であればYをXとして受理

f:id:aidiary:20140712203234p:plain
f:id:aidiary:20140712203238p:plain

ここで、U*M*g(Y) は、[0, M*g(Y)]までの一様乱数を意味している。これが、f(Y) 以下であれば採用するので目標分布の下側にサンプルがきた場合に提案分布fの乱数として Y が採用されることになる。

なるほど！横方向にランダムに位置 Y を決めて、次に縦方向にランダムに位置を決め目標分布の下側に位置すればそのサンプルYを採用していたのか。このようにすると目標分布の山の高いところほど採用されるサンプルが多くなり、山の低いところほど採用されるサンプルが少なくなるはずだ。実際に採用されるサンプルと棄却されるサンプルを色違いで描いてみた。

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

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

# 提案分布g
# 提案分布から乱数生成するためにgvも保持
gv = uniform
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, 3.0))

plt.show()