こんにちは。mutopsyです。この記事では,StanやJAGSのfitオブジェクトを渡すとパラメータのMAP推定値とHDI (最高密度区間) をまとめて返してくれる自作関数を紹介します。
StanやJAGSでパラメータの推定結果の要約を出力すると,通常は事後分布のEAP (i.e., 平均) や分位点の値が出力されますが,事後分布のMAP推定値やHDIが知りたい場合もありますよね。そんなわけで,そういう関数を作ってみました (もう既にあるかもしれないけど知ったこっちゃない)。自作関数といっても,HDIの計算にはHDIntervalというパッケージを利用しますが。関数は以下の通り。
Rコード01: MCMCサンプルからMAP推定値とHDIを計算する関数
library(HDInterval)
maphdi <- function(fit,credMass=0.99,include=TRUE,pars=NA){
if(credMass<0|credMass>1) print("警告:credMassの値が正しくありません。")
Class <- class(fit)[1]
if(Class=="stanfit"){
e <- as.data.frame(rstan::extract(fit))
}else if(Class=="rjags"){
e <- as.data.frame(fit$BUGSoutput$sims.list)
}else if(Class=="data.frame"){
e <- fit
}else{
print("エラー:未対応なオブジェクトです。")
stop()
}
if(is.na(pars[1])) pars <- colnames(e)
if (include==TRUE){
e <- e[,pars]
}else{
e <- e[,!colnames(e) %in% pars]
}
if(length(pars)==1){
MAP <- density(e)$x[which.max(density(e)$y)]
HDI95 <- hdi(e)
HDIAny <- hdi(e,credMass = credMass)
Result <- data.frame(MAP = MAP, Lower95 = HDI95[1], Upper95 = HDI95[2],
LowerAny = HDIAny[1],UpperAny = HDIAny[2])
rownames(Result) <- pars
}else{
MAP <- apply(e,2,function(z) density(z)$x[which.max(density(z)$y)])
HDI95 <- apply(e,2,HDInterval::hdi)
HDIAny <- apply(e,2,function(x) HDInterval::hdi(x,credMass = credMass))
Result <- data.frame(MAP = MAP, Lower95 = HDI95[1,], Upper95 = HDI95[2,],
LowerAny = HDIAny[1,],UpperAny = HDIAny[2,])
}
colnames(Result)[4:5] <- c(
paste("Lower",round(credMass*100,2),sep=""),
paste("Upper",round(credMass*100,2),sep="")
)
return(Result)
}
試しにRからStanのfitオブジェクトを渡してみます。データとモデルは何でもよいので適当なものを渡します。
Rコード02: 関数を使ってみる
>fit.eg <- sampling(stanmodel, data=data, seed=123,
chains=4, iter=10500, warmup=500, thin=1)
> maphdi(fit.eg)
MAP Lower95 Upper95 Lower99 Upper99
beta.1 1.8502036 1.3540767 2.3442492 1.1987086 2.5187496
beta.2 0.6478270 0.5188389 0.7805164 0.4767612 0.8295815
beta.3 0.7103451 0.5996639 0.8201469 0.5627496 0.8524817
beta.4 -0.5622637 -0.7994808 -0.3030058 -0.8891769 -0.2392025
sigma 0.3135227 0.2814639 0.3548575 0.2728315 0.3686322
lp__ 97.7211912 93.6796922 99.1662181 91.7998798 99.3115815
こんな風に,Stanのfitオブジェクトをmaphdi()に引数として渡すだけで各パラメータのMAP推定値と95% HDI・99% HDIを計算して出力することができます。Stan2.18, rstan2.18.2, Rtools3.5, R3.5.2, Windows10で動作確認済み。JAGSにも対応しています (JAGS4.3.0, R2jags0.5-7で確認済み)。
また,fitオブジェクトの代わりに,MCMCサンプルが格納されたデータフレームを渡してもMAP推定値とHDIを返してくれます。
Rコード03: MCMCサンプルを直接渡す
> e <- as.data.frame(rstan::extract(fit.eg))
> head(e)
beta.1 beta.2 beta.3 beta.4 sigma lp__
1 2.194799 0.5785987 0.6255656 -0.4086454 0.3094092 97.46736
2 1.726557 0.6609587 0.7787883 -0.7231942 0.3058624 97.31513
3 2.053161 0.6021407 0.6534262 -0.3992853 0.3425830 96.40984
4 1.941408 0.6503287 0.6399921 -0.4178822 0.2976153 97.93040
5 2.518750 0.4448252 0.7676119 -0.7428853 0.3209268 91.90027
6 1.758346 0.6890341 0.6522970 -0.4046868 0.2985565 97.78386
> maphdi(e)
MAP Lower95 Upper95 Lower99 Upper99
beta.1 1.8502036 1.3540767 2.3442492 1.1987086 2.5187496
beta.2 0.6478270 0.5188389 0.7805164 0.4767612 0.8295815
beta.3 0.7103451 0.5996639 0.8201469 0.5627496 0.8524817
beta.4 -0.5622637 -0.7994808 -0.3030058 -0.8891769 -0.2392025
sigma 0.3135227 0.2814639 0.3548575 0.2728315 0.3686322
lp__ 97.7211912 93.6796922 99.1662181 91.7998798 99.3115815
引数のcredMassに0〜1の範囲の値を与える事で,99% HDIの代わりに任意の確度のHDIを出力することもできます。95% HDIは必ず出力されるようになっています。
Rコード04: 任意のHDIを出力
> maphdi(e,credMass = 0.50)
MAP Lower95 Upper95 Lower50 Upper50
beta.1 1.8502036 1.3540767 2.3442492 1.7011924 2.0385575
beta.2 0.6478270 0.5188389 0.7805164 0.6060299 0.6964437
beta.3 0.7103451 0.5996639 0.8201469 0.6729364 0.7490473
beta.4 -0.5622637 -0.7994808 -0.3030058 -0.6429855 -0.4713128
sigma 0.3135227 0.2814639 0.3548575 0.3030061 0.3280369
lp__ 97.7211912 93.6796922 99.1662181 96.7857008 98.5912001
引数のparsにパラメータの名前を与えれば,特定のパラメータの結果のみを出力することもできます。また,include = FALSE とすると,指定したパラメータ以外の結果を出力することもできます。
Rコード05: 一部のパラメータのみ計算
> maphdi(fit.eg,pars=c("beta.1","beta.2"))
MAP Lower95 Upper95 Lower99 Upper99
beta.1 1.850204 1.3540767 2.3442492 1.1987086 2.5187496
beta.2 0.647827 0.5188389 0.7805164 0.4767612 0.8295815
> maphdi(fit.eg,pars=c("beta.3"))
MAP Lower95 Upper95 Lower99 Upper99
beta.3 0.7103451 0.5996639 0.8201469 0.5627496 0.8524817
> maphdi(fit.eg,include=FALSE,pars=c("beta.1","beta.2"))
MAP Lower95 Upper95 Lower99 Upper99
beta.3 0.7103451 0.5996639 0.8201469 0.5627496 0.8524817
beta.4 -0.5622637 -0.7994808 -0.3030058 -0.8891769 -0.2392025
sigma 0.3135227 0.2814639 0.3548575 0.2728315 0.3686322
lp__ 97.7211912 93.6796922 99.1662181 91.7998798 99.3115815
以上です。Enjoy!