棄却サンプリング - 機械学習に詳しくなりたいブログ

確率分布 $p(x)$ からサンプリングしたいが、正規分布や一様分布のようにライブラリが提供されていない場合にどうするか？逆変換法もその１つですが、棄却サンプリングはより直感的な手法です。

確率分布は $\displaystyle \int p(x)dx = 1$ を満たすべきだが、正規化定数が不明で、

$\displaystyle p(x) = \frac{1}{Z_p}\tilde{p}(x) \tag{1}$

の $\tilde{p}(x)$ のみ既知である場合を考えます。もちろん正規化定数も既知でも問題ないのですが、おそらく不明の場合のほうが一般的だろうと思われます。

全ての $x$ に対して

$kq(x) \ge \tilde{p}(x) \tag{2}$

を満たす定数 $k$ および提案分布 $q(x)$ を定めます。視覚的には、提案分布が $p(x)$ を覆っていれば良いという感じですね。この提案分布 $q(x)$ からはサンプリングが可能なものを選択します。棄却サンプリングの手順は、

$q(x)$ から $z$ をサンプリングする
$\displaystyle \frac{\tilde{p}(z)}{kq(z)}$ の確率で $z$ を採択し、採択できなければ再び $q(x)$ から新たなサンプルをとり、採択できるまで繰り返す。

です。手順2によって採択できない、つまり棄却されるサンプル数を少なくするため、提案分布は $p(x)$ に近いことが望ましい。また、手順２で確率として解釈できるために、つまり $[0,1$ ]の範囲であるために式(2)の条件が必要なんですね。

非常にシンプルでわかりやすい手順です。ああ、こんなことでいいんだって感じ。

さてコードを書いて実験してみました。 $p(x) =\sin^{2} 0.5 \pi x$ からサンプリングしてみます。区間は $[-1,1$ ]とします。わかりやすく提案分布は一様分布としました。つまり $q(x)=0.5$ です。 $p(x)$ の最大値は1になるので、定数 $k=2$ とします。以下のグラフが結果です。

f:id:opabinia2:20180324224014p:plain

$p(x)$ からうまくサンプリングできていますね。

今回のコードです。なお今回は $p(x)$ が対象区間で積分すると1になるようにしましたが、そうでない場合は、このコードのままグラフを描くとサンプリングが失敗しているように見えます。最初はコードにミスがあると思って悩んでしまいました。下のグラフは $p(x) =\sin^{2} 0.3 \pi x$ とした場合です。 $p(x)$ とサンプルのヒストグラムを1つのグラフに描くと形状一致していないように見えてしまいます。 f:id:opabinia2:20180324224945p:plain

# 棄却サンプリングの実験

import numpy as np
import matplotlib.pyplot as plt


def p(x):
    return np.sin(0.5*np.pi*x)**2


# 棄却サンプリング
def sampling():
    k = 2
    # 採択するまでループ
    while True:
        # 提案分布q(z)からサンプリング
        z = np.random.uniform(-1, 1)
        # [0,kq(z)]の一様分布からサンプリング
        u = k*np.random.uniform(0, 0.5)
        # 棄却するか判定
        if p(z) > u:
            return z


x = np.linspace(-1, 1, 500)
y = p(x)

# サンプル数
N = 100000

samples = np.array([sampling() for i in range(N)])

plt.xlim(-1, 1)
plt.hist(samples, bins=50,normed=True,alpha=0.5)
plt.plot(x, y)
plt.show()