Rのplotを使ってNealの漏斗を可視化してみる

今日は、「StanとRでベイズ統計モデリング」という本の10章にあるNealの漏斗問題の可視化を簡単なコードだけで行ってみたいと思います。

ちなみにStanとRでベイズ統計モデリングとは、こちらの本です。

こちらの本、全てのコードがしっかりサポートページに掲載されていて、とても勉強になるのですが、10章にあるNealの漏斗の図は省略されていて(私が見た限り載っていなくて/2019年1月3日現在)、どうやって可視化すればいいんだろうと思ったのでやってみました。
とりあえず、綺麗な図というよりは確認出来ればいいので簡単なplotで試してみたいと思います。

スポンサーリンク

Nealの漏斗問題

まず、Nealの漏斗問題とは何かというと、ベイズモデルを作った際に、パラメータ同士対数事後周辺尤度が次のような漏斗状の形になることで生じるサンプリング問題を指します。

このような状況になると、愚直にStanでモデリングを行っても、サンプルを効率よく取ってくることが出来ません。そこで再パラメータ化という、一旦独立な分布からサンプリングしてやったのち、変数変換を行う処置をとるわけです。

Rコード

では、コードを打っていきます。

まず、rstanの読み込み

library(rstan)

次にモデルを作っておきます

model="
parameters{
	real a;
	vector[1000] r;
}

model{
	a~normal(0,3);
	r~normal(0,exp(a/2));
}"

fit_a=stan(model_code=model,seed=1234,chains=1,iter=4000)
a1=rstan::extract(fit_a)$a
r1=rstan::extract(fit_a)$r[,1]

この状態でのサンプリングはあまりうまくいきません。それをplotでチェックします。

r=seq(-10,10,0.025)
a=seq(-10,10,0.025)
n=length(r)
lp1=matrix(,n,n)
for(i in r){
	for(j in a){
		lp1[which(r==i),which(a==j)]=log(dnorm(i,0,3))+log(dnorm(j,0,exp(i/2)))
	}
}
lp1[lp1==-Inf]=0
contour(r,a,t(lp1))
points(r1,a1,col="red")

あまりうまくサンプリング出来ていないことが確認出来ました。

さて、再パラメータ化を行った場合はどうでしょう。

model2="
parameters{
	real a_raw;
	vector[1000] r_raw;
}

transformed parameters{
	real a;
	vector[1000] r;
	a=3.0*a_raw;
	r=exp(a/2)*r_raw;
}

model{
	a_raw~normal(0,1);
	r_raw~normal(0,1);
}"

fit_b=stan(model_code=model2,seed=1234,chains=1,iter=4000)
a2=rstan::extract(fit_b)$a
r2=rstan::extract(fit_b)$r[,1]

重ねてみる

contour(r,a,t(lp1))
points(r2,a2,col="green")

漏斗の先も取ってきて来れていることが確認出来ました！