俺なりの事後分布からのサンプリング

前提

ベイズの定理
\[ p(X|Y) = \frac{p(Y|X)p(X)}{p(Y)} \]
を使って、事後分布からのサンプリング、つまりyが与えられている条件付きの下でxからのサンプリングを行いたい。
ここでは確率変数\( X,Y \)が相関付きの2変量（標準）正規分布に従っているとする。

従って、条件付きの確率密度関数はすぐに計算できて
\[ p(x|Y=y)=\frac{1}{\sqrt{2\pi(1-\rho^2)}} e^{-\frac{(x -\rho y)^2}{2(1-\rho^2)}} \]
つまり平均\( \rho y \)・標準偏差\( \sqrt{1-\rho^2} \)の正規分布に従う。

このドキュメントはこの事後分布からのサンプリングを俺なりのやり方でやってみたというお話。

普通にサンプリング

ほとんど自明な話だが、俺なりのやり方の答え合わせ用途も兼ねて。

今回のケースだと解析的に事後分布が計算出来ているのでその分布に従う乱数を生成すれば、\( y \)で条件付けられた下での\( x \)の標本を取得出来る。

\( y \)を1だと条件付し、相関係数\( \rho \)を0.5としてシミュレーションしてみると

rho <- 0.5
y.fixed <- 1
N <- 10000
x <- rnorm(N, rho * y.fixed, sqrt(1 - rho^2))
x.mean <- mean(x)  # 0.5*1 = 0.5
x.sd <- sd(x)  #sqrt(1-0.5^2) = 0.866
c(x.mean, x.sd)

## [1] 0.4822 0.8616

# PLOT
h <- hist(x, sqrt(N), freq = FALSE)
xfit <- seq(min(x), max(x), length = 40)
lines(xfit, dnorm(xfit, rho * y.fixed, sqrt(1 - rho^2)), col = "blue", 
    lwd = 2)

plot of chunk unnamed-chunk-1

となって、大体いい感じにサンプリング出来ている（青線は解析的な事後分布の密度関数グラフ）。

俺なりのサンプリング

今回は２変量正規分布での話なので、上述のように難なく事後分布からのサンプリングが出来たが、実際のデータやもっと別の分布を扱おうと思った場合、こう簡単にはできないケースが多く、実際にはギブスサンプリング等の手法が用いられるわけです。

もちろんそれでもいいんだけど、俺なりのやり方って奴を考えてみたのでここで試しにやってみる。
ベイズの定理
\[ p(X|Y) = \frac{p(Y|X)p(X)}{p(Y)} \]
の右辺をじっと睨んでやって、\( \frac{p(Y|X)}{p(Y)} \)と\( p(X) \)の２つの項に分けて以下のように考えてみた。

俺なりの事後分布からのサンプリングアルゴリズム

\( p(X) \)（今回の場合は標準正規分布）からのサンプルを引く
サンプリングした各\( x_i \)に対して確率密度の比\( \frac{p(Y|X)}{p(Y)} \)を計算しウェイトとする（尤度をウェイトとしてるという意味に近い）
ウェイトを整数にしたいので、適当に切り上げ・スケーリングする。なので2の操作での分母\( p(Y) \)は全体共通なので実際にはいらない
サンプリングした各\( x_i \)を計算したウェイト個倍し、標本集団\( s \)とする。
標本集団\( s \)からブートストラップサンプリングし所望の条件付き分布からのサンプリングとする

意味合いとしては\( p(X) \)から引いたサンプルに条件付きの効果を乗せるために尤度をウェイトとし、リサンプリングしているという事。このやり方だと尤度関数（今回の場合は\( P(Y|X) \)も正規分布だが・・・）だけわかってれば計算できるのでいいんじゃないかと思ってる。
実際上記の手順でサンプリングし、平均・標準偏差・分布形をチェックしてみると

rho <- 0.5
y.fixed <- 1
N <- 10000
x <- rnorm(N)  #事前分布p(X)からのサンプリング
weight <- sapply(x, function(z) dnorm(y.fixed, rho * z, sqrt(1 - 
    rho^2))/dnorm(y.fixed))  #尤度的な物P(Y|X)/P(Y)の計算。分母のP(Y)はなくてもOK
weight <- round(10 * weight)  #スケーリング&整数化
s <- rep(x, weight)  #ウェイトに応じてサンプルを増やす
x <- sample(s, N)  #標本集団sからサンプリングし、条件付分布からのサンプリングとする
x.mean <- mean(x)  #0.5
x.sd <- sd(x)  #0.866
c(x.mean, x.sd)

## [1] 0.4900 0.8585

# PLOT
h <- hist(x, sqrt(N), freq = FALSE)
xfit <- seq(min(x), max(x), length = 40)
lines(xfit, dnorm(xfit, rho * y.fixed, sqrt(1 - rho^2)), col = "blue", 
    lwd = 2)

plot of chunk unnamed-chunk-2

・・・たまに歯抜けデータな所があるのはなんとかしたいが、大体うまく出来ているように見えるよね？（青線は解析的な事後分布の密度関数グラフ）

既に誰かがやってそうなのでこのやり方の手法名とかこの手の話に詳しい書籍等ご存じな親切な方がいらっしゃいましたらTwitter: @teramonagiまで連絡下さい。