ベイズマックスへようこそ。
このブログでは,心理学と統計が大好きなmutopsyが,
統計に関するあれやこれやを自由気ままに投稿します。

2019年01月04日

StanやJAGSの計算結果からパラメータのMAP推定値とHDIをまとめて得るためのRの自作関数

 こんにちは。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!

【Rに関するTipsの最新記事】
posted by mutopsy at 19:10 | Rに関するTips