事後最頻値をとる関数(get_posterior_map)を作りました。

どうも。

この記事は、StanAdventも何も関係ありませんｗ

Stanで事後平均値(EAP)をとってくるときにget_posterior_mean()があると思いますが、事後分布の最頻値(MAP)をとってきたい場合もたまにあるかと思います。

僕はその都度R codeを書いていたのですが、めんどくさくなったので、MAP推定値をとってくる関数を書きました。
それ以上でもそれ以下でもない記事ですｗ

依存パッケージは、dplyrだけです(%>%が使いたいだけ)。

構造：get_posterior_map(stanfit, pars)
stanfit : Stan結果オブジェクト
pars : 取りたいパラメータ

注意：抽出したいパラメータの次元が4次元まで対応しています。

使用したい場合、以下のコードをコピーして実行してください。

get_posterior_map <- function(stanfit,pars) {
  if(!require("dplyr",character.only = T)) install.packages("dplyr", dependencies = T); require("dplyr",character.only = T)
  maxpost <- function(par){ density(par)$x[which.max(density(par)$y)]}
  r <- rstan::extract(stanfit,pars=pars)
  Dim <- r[[pars]] %>% dim %>% length - 1
  if (Dim == 1 ){
    map_est <- array(NA,dim=c(dim(r[[pars]])[2]))
    for (l in 1:dim(r[[pars]])[2])
      map_est[l] <- r[[pars]][,l] %>% maxpost(.)
    return(map_est)
  }else if (Dim == 2 ){
    map_est <- array(NA,dim=c(dim(r[[pars]])[2],dim(r[[pars]])[3]))
    for (l in 1:dim(r[[pars]])[2]) 
      for (m in 1:dim(r[[pars]])[3]) 
        map_est[l,m] <- r[[pars]][,l,m] %>% maxpost(.)
    return(map_est)
  }else if (Dim == 3 ){
    map_est <- array(NA,dim=c(dim(r[[pars]])[2],dim(r[[pars]])[3],dim(r[[pars]])[4]))
    for (l in 1:dim(r[[pars]])[2]) 
      for (m in 1:dim(r[[pars]])[3]) 
        for (n in 1:dim(r[[pars]])[4]) 
          map_est[l,m,n] <- r[[pars]][,l,m,n] %>% maxpost(.)
    return(map_est)
  }else if (Dim == 4){
    map_est <- array(NA,dim=c(dim(r[[pars]])[2],dim(r[[pars]])[3],dim(r[[pars]])[4],dim(r[[pars]])[5]))
    for (l in 1:dim(r[[pars]])[2]) 
      for (m in 1:dim(r[[pars]])[3]) 
        for (n in 1:dim(r[[pars]])[4]) 
          for (o in 1:dim(r[[pars]])[5]) 
            map_est[l,m,n,o] <- r[[pars]][,l,m,n,o] %>% maxpost(.)
    return(map_est)
  }else{
    print("Check Parameter Dimensions")
  }
}

って感じです。
もっとうまく書く方法はあると思いますが、僕はこれでうまく取ってこれるので問題無いです:)

Enjoy！