ベイズ統計学 | 事後分布のサンプリング Rejection Sampling

サンプリングが難しい事後分布を近似的にサンプリングするための方法を紹介します．今回はRejction Samplingです．この方法はMCMCの基礎となるAcceptance-Rejection Methodの1つです．

前回は事後分布を決定的な方法で近似するDirect Methodを紹介しました．

サンプリングしたい事後分布$p(\theta | y)$のことを目標分布 (target distribution) と呼びます．

Rejection Samplingのアルゴリズム

さて，Rejection Samplingの簡単な概略を示します．まずサンプリングが簡単にできる$g(\theta)$を考えます．

ただし，$g(\theta)$は目標分布$p(\theta|y)$のサポート（取りうる$\theta$の範囲）より大きい必要があります．それだけでなく，importance ratioである$p(\theta | y)/g(\theta)$が上に有界である必要があります．つまり，ある上界 (upper bound) $M$が存在して

$$
\frac{p(\theta | y)}{g(\theta} \leq M.
$$

実はこの条件は，$g(\theta)$の方がサポートが大きい（か同じ）という条件を含むので上界$M$が存在するような$g(\theta)$を見つければ良いです．

さて$g(\theta)$を決定した後は，$g(\theta)$から$\theta^s$をサンプリングします．そして，そのサンプル$\theta^s$を$p(\theta^s|y)/Mg(\theta^s)$の確率でAcceptします（$1 - p(\theta^s|y)/Mg(\theta^s)$の確率でRejectします）．

Accept-Rejectの判断は$[0, 1]$の一様分布の乱数$U$を発生させて，$U <p(\theta^s|y)/Mg(\theta^s)$であればAcceptし，そうでなければRejectします．

$\theta^s$は$g(\theta)$からサンプルされます．全てAcceptしてしまえば$\theta^s$たちは$g(\theta)$からのサンプルです．しかし，$p(\theta|y)/Mg(\theta)$の確率でAcceptするという条件付けを行うことで，$p(\theta | y)$からのサンプルになります．

以下の絵を見るとなんとなく納得するかもしれません．

次に証明を与えます．

Rejection Samplingの証明

$p(\theta | y)$，$g(\theta)$の累積分布関数をそれぞれ$F(\theta), G(\theta)$とおきます（観測データ$y$の条件付けは省略します）．また$X, Z$をそれぞれ$F, G$に従う確率変数とします．

$[0, 1]$の一様分布の確率変数を$U$とおきます．示すべきはAcceptされたという条件付の$Z$の確率変数が$X$と等しくなることです．つまり，

$$
F(\theta) = \Pr \big(Z \leq \theta | U < \frac{p(Z | y)}{Mg(Z)}\big).
$$

まずAcceptされる確率$\Pr(U \leq p(Z|y)/Mg(Z))$を計算します．

$$
\begin{aligned}
\Pr(U \leq p(Z|y)/Mg(Z)) & =  \int_{-\infty}^{\infty} \int_0^1 \frac{p(z|y)}{Mg(z)} g(z) dude \\
& = \frac1M
\end{aligned}
$$

次に$\Pr \big(U < \frac{p(Z | y)}{Mg(Z)}| Z \leq \theta \big)$を計算します．

$$
\begin{aligned}
\Pr \big(U < \frac{p(Z | y)}{Mg(Z)}| Z \leq \theta \big) & = \frac{\Pr \big(U < \frac{p(Z | y)}{Mg(Z)}, Z \leq \theta \big)}{G(\theta)} \\
& = \frac1{G(\theta)} \int_{-\infty}^\theta  \Pr \big(U < \frac{p(w | y)}{Mg(w)} | Z = w \big)g(w)dw \\
& =  \frac1{G(\theta)} \int_{-\infty}^\theta   \frac{p(w | y)}{Mg(w)} g(w)dw \\
& = \frac{F(\theta)}{MG(\theta)}
\end{aligned}
$$

上で計算したものを統合すると，

$$
\begin{aligned}
\Pr \big(Z \leq \theta | U < \frac{p(Z | y)}{Mg(Z)}\big) & = \frac{\Pr \big(U < \frac{p(Z | y)}{Mg(Z)}| Z \leq \theta \big)G(\theta)}{\Pr(U < \frac{p(Z | y)}{Mg(Z)})} \\
&= F(\theta).
\end{aligned}
$$

Rでシミュレーション

Rejection Samplingを実際にRで実装してみます．

今回は二項分布を尤度 (likelihood) としたモデルを考えます．

$$
y | p \sim \text{Bin}(n, p)．
$$

サンプルサイズ$n$は$20$とします．また真の成功確率のパラメタ$p^*$は$0.4$であるとします．真のパラメタを持った二項分布から観測データである$y$をサンプルします．

また$p$の事前分布（prior distribution) はベータ分布とします．そのハイパーパラメータは$\alpha = \beta = 1/2$とします．

このとき，$p$の事後分布 (posterior distribution) は$\text{Be}(y+\alpha, n - y + \beta)$です．実際はベータ分布はRの関数を使えば簡単に（疑似）乱数を生成できます．

次に$g(p)$を決めなければいけません．今回は$[0, 1]$の一様分布とします．次にMを決めますが，$[0, 1]$の一様分布の場合は$g(p) = 1$なので事後分布の最大値を$M$として使えます．つまり，$\text{Be}(y+\alpha, n - y + \beta)$の最頻値 (mode) における密度関数の値を使えば良いので，$p_{\text{max}} = \frac{y + \alpha -1}{n + \alpha + \beta - 2}$で評価した密度関数の値がMとなります．

実際のコードが以下のようになります．

n <- 20
a <- b <- 0.5
p <- 0.4
y <- rbinom(1, size=n, prob = p)
iter <- 100000

post.a <- a + y
post.b <- n - y + b
p.max <- (post.a - 1) / (post.a + post.b - 2)
M <- dbeta(p.max, post.a, post.b)

p.sims <- 0
j <- 0
for (i in 1:iter) {
 p.cand <- runif(1)
 r <- dbeta(p.cand, post.a, post.b) / (M * 1)
 u <- runif(1)
 if (u < r) { # accept
   p.sims[j] <- p.cand
   j <- j + 1
 }
}
hist(p.sims ,breaks=500, xlab='p', xlim=c(0,1), main='Sample obtained using rejection sampling', freq=FALSE)
curve(dbeta(x, post.a, post.b),xlab='p',ylab='p(p|y)',add=TRUE,col=2,lwd=2)

Rejection Samplingの弱点

都合の良い$g(\theta)$を見つける必要があります．上で考えたような二項分布モデルでは$n$が大きくなるにつれて，最大値が大きくなり，いつまでも一様分布をサンプルの分布として用いるとAcceptされるサンプルサイズが小さくなり非効率化してしまいます．漸近理論に基づいてtruncatedされた正規分布を用いることもあるようです．

Reference

Gelman, Andrew, John B. Carlin, Hal S. Stern, and Donald B. Rubin. Bayesian data analysis. Chapman and Hall/CRC, 1995.