Metropolis-Hastings（MH）アルゴリズムを用いた正規母集団パラメータのベイズ推定

Xi （i = 1～n）を正規母集団N（μ、σ²）からのi.i.d.サンプルとする。MHアルゴリズムによるμとσ²のベイズ推定例を示す。
（参考文献）大森裕浩「マルコフ連鎖モンテカルロ法の基礎と統計科学への応用」、計算統計II（第III部）、岩波書店、p.168。

１．事前分布

２．尤度関数

３．事後分布

４．未知パラメータの全条件付事後分布

５．MHアルゴリズム
各パラメータの全条件付事後分布に基づくアルゴリズムを採用する。提案分布が対称（q(a,b)=q(b,a)）のとき、MHアルゴリズムはMetropolisアルゴリズムと呼ばれる。また、提案分布が全条件付事後分布に等しい場合は、下の採択確率が1となり、p(μ'|σ_(k)²,x）及びp((σ²)'|μ_(k+1),x）から順次サンプリングを行うことになるので、MHアルゴリズムはGibbsサンプラーに一致する。

６．Rコードの一例
#μとσ²ともに正規分布を提案分布として用いた場合の例。
library(MCMCpack)
X <- rnorm(100,5,1)
n <- NROW(X)
Xbar <- mean(X)
Sum2 <- sum((X-Xbar)^2)
mu0 <- 0.0 # ハイパーパラメータ
sig20 <- 1000 # ハイパーパラメータ
nu0 <- 0.0005 # ハイパーパラメータ
ramda0 <- 0.0005 # ハイパーパラメータ
N <- 100000 # サンプル数
Nb <- 1000 # burn-in
sample_mu <- 1:N
sample_sd <- 1:N
sample_var <- 1:N
sig2 <- 1.0 # 初期値
sig <- 1.0 # 初期値
mu <- 0 # 初期値
var <- 1 # 初期値
N_mu <- 0 # 採択率計算用
N_var <- 0 # 採択率計算用
for(i in -Nb:N) {
mean <- (n*Xbar/sig2+mu0/sig20)/(n/sig2+1/sig20) # 直前のsig2を用いる。
sd <- sqrt(1/(n/sig2+1/sig20))
a1 <- Xbar # 提案分布のパラメータ ： 独立連鎖
b1 <- 0.5 # 提案分布のパラメータ（採択率の調整）
cand_mu <- rnorm(1,a1,b1)
alpha <- min(1, (dnorm(cand_mu, mean, sd)/dnorm(mu, mean, sd))*
(dnorm(mu,a1,b1)/dnorm(cand_mu,a1,b1)))
u <- runif(1)
if(u < alpha){
mu <- cand_mu
if(i > 0) N_mu <- N_mu+1 #採択数カウント
}
nu1 <- (n+nu0)/2
ramda1 <- (Sum2+n*(mu-Xbar)^2+ramda0)/2 # 直前のmuを用いる。
a2 <- 1 # 提案分布のパラメータ ： 独立連鎖
b2 <- 0.1 # 提案分布のパラメータ（採択率の調整）
cand_var <- rnorm(1,a2,b2)
alpha <- min(1, (dinvgamma(cand_var,nu1,ramda1)/dinvgamma(var,nu1,ramda1))*
(dnorm(var,a2,b2)/dnorm(cand_var,a2,b2)))
u <- runif(1)
if(u < alpha){
var <- cand_var
if(i > 0) N_var <- N_var+1 # 採択数カウント
}
sig2 <- var # 初めに戻ってmeanの計算に用いる。
sig <- sqrt(var)
if(i > 0) {
sample_mu[i] <- mu #μのサンプル
sample_var[i] <- var #σ^2のサンプル
sample_sd[i] <- sig #σのサンプル
}
}
N_mu/N # 採択率
N_var/N # 採択率
summary(sample_mu)
summary(sample_var)
st_mu <- c(mean(sample_mu),sd(sample_mu))
st_var <- c(mean(sample_var),sd(sample_var))
print(st_mu)
print(st_var)
plot(sample_mu,type="l")
par(ask=TRUE)
plot(density(sample_mu))
par(ask=TRUE)
plot(sample_var,type="l")
par(ask=TRUE)
plot(density(sample_var))

７．結果の表示
> N_mu/N #採択率
[1] 0.25258
> N_var/N #採択率
[1] 0.89771
> summary(sample_mu)
Min. 1st Qu. Median Mean 3rd Qu. Max.
4.623 5.027 5.095 5.095 5.163 5.485
> summary(sample_var)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.5566 0.9185 1.0090 1.0220 1.1120 1.6640
> st_mu <- c(mean(sample_mu),sd(sample_mu))
> st_var <- c(mean(sample_var),sd(sample_var))
> print(st_mu)
[1] 5.0947409 0.1003362
> print(st_var)
[1] 1.0220987 0.1463964