事後最頻値をとる関数(get_posterior_map)を作りました。
どうも。
この記事は、StanAdventも何も関係ありませんw
Stanで事後平均値(EAP)をとってくるときにget_posterior_mean()があると思いますが、事後分布の最頻値(MAP)をとってきたい場合もたまにあるかと思います。
僕はその都度R codeを書いていたのですが、めんどくさくなったので、MAP推定値をとってくる関数を書きました。
それ以上でもそれ以下でもない記事ですw
依存パッケージは、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!