マルコフ連鎖モンテカルロ法（MCMC）には様々なアルゴリズムがあり、メトロポリス・ヘイスティングス法は尤度関数があれば推定ができる使い勝手のよい方法だと知られています。しかし、高次元で効率が悪いことが知られており、また多峰の確率分布には上手く対応できないそうで、最近のベイズ推定ではHamiltonian Monte Carlo法や、その発展であるNo-U-Turn Samplerが使われています。
Hamiltonian Monte Carlo法は教科書的な手法のようで、Rにもhmclearnと言うHamiltonian Monte Carlo法の教育用パッケージがあるのですが、参考にしつつFortranで書いてみました。推定モデルにLinear Mixed EffectsがなくOrdered Logitがあったり、初期値とステップ幅を最尤法の結果から設定するようにするなど*1、細部では変更があります。
https://github.com/uncorrelated/HMCFortran
Linux環境でgit cloneしたディレクトリでR CMD INSTALL . でRにインストールして使うようになっていますが、Fortranのライブラリを2つ（と片方に内蔵されたCのライブラリを1つ）をリンクします。つまり、lbfgsbとdSFMT_F03_interfaceを、どこかにディレクトリー（e.g. /var/tmp）を作ってgit cloneして、MakefileのSHSRC=/var/tmpを適時書き換える必要があります*2。
cd /var/tmp git clone https://github.com/uncorrelated/HMCFortran.git git clone https://github.com/jacobwilliams/lbfgsb.git git clone https://github.com/jsspencer/dSFMT_F03_interface.git cd HMCFortran R CMD INSTALL .
使い方は以下のように、なるべくR的にしました。モデル式を書いて、summaryやplotなどのジェネリック関数で結果が見られるようにしています。
library(hmcfortran) # data generator for practice mkdata <- function(nr, nc, sd = 1){ X <- matrix(c(rep(1, nr), runif(nr*(nc-1), -2, 2)), nr) colnames(X) <- c("Const.", letters[1:(nc - 1)]) beta <- c(0.5, rep(1, nc - 1))*(-1)^(1:nc) y <- X %*% beta + rnorm(nr, sd = sd) cbind(y, as.data.frame(X)) } np <- 4 # number of parameters data <- mkdata(150, np - 1) # prior/hyper parameters ig.alpha <- 1 ig.gamma <- 1 beta.mean <- rep(0, np - 1) beta.Sigma <- diag(1e+10, np - 1, np - 1) # do Hamiltonian Monte Carlo method to estimate a liner model r_hmc_lm <- hmc.blm(y ~ a + b, ig.alpha, ig.gamma, beta.mean, beta.Sigma, data = data, N = 3000, nchains = 3) # get the marginal likelihood mlf_1 <- log_ml(r_hmc_lm) # alternative model: third parameter, a, is zero. r_hmc_lm_0 <- hmc.blm(y ~ b, 1, 1, rep(0, np - 2), diag(1e+10, np - 2, np - 2), data = data, N = 3000) # get the marginal likelihood of alternative model mlf_0 <- log_ml(r_hmc_lm_0) cat("Laplace Approximated BF_{10}:", exp(mlf_1 - mlf_0), "\n") BF <- exp(log_BF(r_hmc_lm, 3, 0)) cat("IWAMDE Savage-Dickey ratio BF_{01}:", BF, " BF_{10}:", 1/BF, "\n") # show the prediction, where explanatory variables are the mean of X. prdct <- predict(r_hmc_lm, X = apply(r_hmc_lm$input$X, 2, mean)) hist(prdct, breaks = 20) # call coda::plot plot(r_hmc_lm) # calcurate R-hat gelman.diag(r_hmc_lm) # call other functions of coda package # gelman.plot(r_hmc_lm$mcmc.list) # geweke.plot(r_hmc_lm$mcmc.list) # crosscorr.plot(r_hmc_lm$mcmc.list) # autocorr.diag(r_hmc_lm$mcmc.list) # autocorr.plot(r_hmc_lm$mcmc.list) # call coda::summary summary(r_hmc_lm) # call coda::HPDinterval HPDinterval(r_hmc_lm) # statistics of the sample of MCMC mean(r_hmc_lm) median(r_hmc_lm) quantile(r_hmc_lm, 0.5)
線形回帰、二項ロジット、順序ロジット、混合ロジット、ポアソン回帰の他、Rのユーザー定義の目的関数を使えるようにしています。プリセットのモデルでは周辺尤度のラプラス近似と、MCMCサンプルから制約モデルのベイズファクター（IWMDE/Chen (1992) Savage-Dickey density ratio）が計算でき、codaパッケージを利用することでMCMCの収束判定も教科書的に可能なので、MCMCpackと同程度に気軽にベイズ推定ができます。ドキュメントがまったく整っていないので、何とかしたいですね！
今後、気力と体力があればLinear Mixed Effectsと欠損値処理は実装したいのですが、どうも無さそうなので期待しないでください('-' )＼(--；)BAKI
