MCMCpackのMCMCmetrop1RでBayes Factorを雑に計算
マルコフ連鎖モンテカルロ法(MCMC)を使いたいベイジアンの統計ユーザーはStanを使っている人が多く、R使いもフロントエンドRStanを通して利用していることがほとんどだと思うのですが、MCMCpackと言うのもあります。このMCMCpackでベイズ因子(Bayes Factor)を使ってモデル選択をしてみましょう。
MCMCpackにはよく使われるらしい特定の共役事前分布*1を置ける線形回帰やポアソン回帰などの関数が用意されており、それらはベイズ因子を簡単に出力できるように設計されています。しかし、今回は正規分布と対数正規分布のどちらが当てはまりがよいかモデル選択をするので使えません。そして、汎用的なベイズ因子の計算関数は無いようです。
詰まった感があるのですが、推定モデルごとに周辺尤度を計算できればベイズ因子は簡単に計算できます。そして、事後分布の推定が終わっていれば、ベイズの定理から周辺尤度を簡単に計算できます。精度を考えなければ。
ベイズの定理の事後確率と周辺尤度を入れ替えてみましょう。
はデータと言うか観測値のベクトル、はパラメーターのベクトル、は確率、は条件付き確率です。は任意のですが、精度の都合でMAP推定量など、なるべく大きな値を選びます。はで計算した尤度、はの事前確率です。
1. データセット
今回は以下の対数正規分布からの乱数を使います。
set.seed(20211204) n <- 15 y <- rlnorm(n, meanlog=1.2, sdlog=0.8)
真のパラメーターが分かっている方が、推定がうまくいっているか分かりやすいので。
ヒストグラムを描くと、こんな感じです。
2. 対数尤度関数
# 正規分布の対数尤度関数 llf_nd <- function(theta){ if(any(theta < 0.0)){ return(-Inf) } sum(dnorm(y, mean=theta[1], sd=theta[2], log=TRUE)) } # 対数正規分布の対数尤度関数 llf_lnd <- function(theta){ if(any(theta < 0.0)){ return(-Inf) } sum(dlnorm(y, meanlog=theta[1], sdlog=theta[2], log=TRUE)) }
3. MCMCで推定
MCMCpackをインストールしてあったら、簡単にできます。なお、MCMCmetrop1R関数は初期値をoptimでMAP推定してからMCMCを開始するため、最適化アルゴリズムがMAP推定に失敗するとエラーを出します。そしてデフォルト指定のBFGSがエラーを出してきたので、L-BFGS-Bに代えました。
init_theta <- c(mean(y), sd(y)) library(MCMCpack) # 正規分布のモデルの事前分布は(0 7)と(0 5)の一様分布 prior.nd <- function(theta){ dunif(theta[1], min=0, max=7, log=TRUE) + dunif(theta[2], min=0, max=5, log=TRUE) } post.samp.nd <- MCMCmetrop1R(function(theta){ llf_nd(theta) + prior.nd(theta) }, theta.init=log(init_theta), optim.method="L-BFGS-B") # 対数正規分布のモデルの事前分布は(0 log(7))と(0 log(5))の一様分布 prior.lnd <- function(theta){ dunif(theta[1], min=0, max=log(7), log=TRUE) + dunif(theta[2], min=0, max=log(5), log=TRUE) } post.samp.lnd <- MCMCmetrop1R(function(theta){ llf_lnd(theta) + prior.lnd(theta) }, theta.init=log(init_theta))
4. 中央値から周辺尤度を計算
今回は手軽に中央値を使います。ノンパラメトリックな確率密度関数の推定法…とかっこよく言いたいところですが、MCMCを繰り返して精度をあげるChib (1995)のような事はしていないので雑です。
# 中央値を得る medians.nd <- summary(post.samp.nd)$quantiles[, "50%"] medians.lnd <- summary(post.samp.lnd)$quantiles[, "50%"] # パラメーターyに対応する事後確率を求める get_posterior_probability <- function(post.samp, y, nos=100){ p <- length(y) # 次元 nr <- length(post.samp[, 1]) count <- function(a){ v <- 1 flag <- rep(TRUE, nr) for(column in 1:p){ s <- post.samp[, column] range <- max(s) - min(s) # y±range*aに含まれる点を選ぶ flag <- flag & s <= y[column] + range*a & s >= y[column] - range*a v <- v * (2*range*a) } c("nflag"=sum(flag), "volume"=v) } # 区間に入る点がnos個になるように、区間を調整する l <- 0 h <- 1 m <- 0.5 nz <- h for(i in 1:64){ m <- (l + h)/2 n <- count(m)[["nflag"]] - nos if(n == 0){ break; } else if(n < 0){ l <- m } else { h <- m nz <- m } } r <- count(ifelse(0<=n, m, nz)) # r["nflag"]/length(post.samp)は確率(累積分布)になるので、区間の大きさr["volume"]で割って(平均)確率密度にする r[["nflag"]]/length(post.samp)/r[["volume"]] } # 中央値の事後確率 posterior_p.nd <- get_posterior_probability(post.samp.nd, medians.nd) posterior_p.lnd <- get_posterior_probability(post.samp.lnd, medians.lnd) # 対数周辺尤度を作る mll.nd <- llf_nd(medians.nd)-log(posterior_p.nd) + prior.nd(medians.nd) mll.lnd <- llf_lnd(medians.lnd)-log(posterior_p.lnd) + prior.lnd(medians.lnd)
5. ベイズ因子を計算
割り算と言うか、引き算して指数をとって終わりです。
sprintf("Bayes Factor: %.3f", exp(mll.lnd-mll.nd))
8.967でSubstantialと言えるような値になりますが、ヒストグラムの見た目にあった結果ではないでしょうか。なお、シャピロ・ウィルク検定をかけると10%有意で正規分布が棄却されて、対数正規分布は棄却されません。
shapiro.test(x=y)$p.value shapiro.test(x=log(y))$p.value
6. 注意事項
ほぼ思いつきで書いた胡散臭い方法なので、実用の参考にはしないでください*2。