イントロ
この記事はベイズ予測分布を応答変数データのヒストグラムに重ねてplotする方法のメモです．
目的
応答変数が離散的なスケールで測定されている場合，bayesplotのppc_dens_overlay関数が使えなくて面倒くさいってことがあります．
無理矢理やると，こんな感じになり，泣けてきます．
というわけで，離散的な応答変数に連続関数の予測分布を重ねる方法を考えました．
bayesplotの設定で対応できるのかもしれませんが，どうせ後で凡例などを修正したくなるはずなので，先を見越してggplotで自作します．
準備
まずStanの出力stanfitオブジェクトからパラメータの事後分布を取り出し，
ベイズ予測分布をRで定義します．
fit <- extract(stanfit) predict<-function(x) { mean(dnorm(x,fit$mu,fit$sigma)) #確率モデルにパラメータをベクトルで代入 }
確率モデルは正規分布を仮定しています．
ここでfit$mu, fit$sigmaの長さはMCMCサンプルの個数と一致します．たとえばチェイン3で反復1500，バーンインが500と仮定すれば
MCMCサンプルの個数は$(1500-500)\times 3=3000$です．
関数predictは，確率モデルに事後分布の実現値（MCMCサンプル）を代入した値の平均なので，実質的には確率モデルの事後分布による平均（ベイズ予測モデル）の近似になっています（予測分布が解析的に出てくる場合はこの関数は必要ありません）．
次に予測分布predictを使って，pdfの値を計算します．
pdf<-map_dbl(seq(0.10,0.02),predict)
seqの範囲はpdfの長さがMCMCではなくデータのサンプルサイズに一致するように適当に調整します．
なぜかといえば，ggplotで応答変数のヒストグラムとpdfの曲線を重ねたいからです．
ggplotのコ－ドは例えば次のようなものになるでしょう．
d1 <- data.frame(Level=y,yp=pdf,xaxis=seq(0,10,0.02)) ggplot(d1,aes(x=Level,y=..density..)) + geom_histogram(binwidth=1,fill="white",colour="black")+ geom_line(aes(x=xaxis,y=pdf))
結果はたとえば次のとおりです．
テスト用にサンプルサイズを小さくしたので，あまりフィットしていませんが，まあいいでしょう．
説明変数があるとき
さていろいろ試しているうちに，説明変数がある場合の予測分布の定義がよく分からないことに気づきました． パラメータが説明変数によって変化するモデル，典型的には次のような線形回帰モデルです*1．
$$ Y_i \sim N( \beta_{0} +\beta_1 X_i, s^{2}) \quad i=1,2,...,n $$
予測分布の定義を次のように書き換えます．
predict<-function(y) { dnormlist <- c() for(i in 1:n){ dnormlist <- append(dnormlist,dnorm(y,fit$b0+fit$b1*x1[i])) } mean(dnormlist) }
ここでx1が説明変数です．サンプル$i$ごとに事後分布で確率モデルを平均し，さらに最後にサンプル方向に平均します．
テスト用に応答変数が説明変数によって二峰型になるパタンを考えてみましょう
n <- 150 # sample size x1 <- c(rep(0,n/2),rep(5,n/2))#説明変数 y <-c() for(i in 1:n){ y <- append(y,rnorm(1,x1[i],1)) }
予測分布を当てはめます．
多分もっと効率のよい書き方があると思いますので気づいたら直します．
