Gibbsサンプラーを用いた正規母集団パラメータのベイズ推定

Xi （i = 1～n）を正規母集団N（μ、σ²）からのi.i.d.サンプルとする。Rを用いたGibbsサンプラーによるμとσ²のベイズ推定例を示す。
（参考文献）中妻照雄著「入門ベイズ統計学」、朝倉書店、p.145。

１．事前分布

２．尤度関数

３．事後分布

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

５．Gibbsサンプラーのアルゴリズム

６．Rコードの一例
library(MCMCpack) # rinvgamma（）を使用するため
X <- c(0.0459, 0.0436, 0.0207, 0.0867, 0.1678, 0.1748) #データ
n <- NROW(X)
Xbar <- mean(X)
Sum2 <- sum((X-Xbar)^2)
mu0 <- 0.0 # ハイパーパラメータ
sig20 <- 0.01 # ハイパーパラメータ （sig20はτ₀²を表す。）
nu0 <- 36.0 # ハイパーパラメータ
ramda0 <- 1.36 # ハイパーパラメータ
N <- 100000 # サンプル数
Nb <- 1000 # burn-in
sample_mu <- 1:N
sample_sd <- 1:N
sample_var <- 1:N
sig2 <- 1.0 # 初期値
for(i in -Nb:N) {
mean <- (n*Xbar/sig2+mu0/sig20)/(n/sig2+1/sig20)
sd <- sqrt(1/(n/sig2+1/sig20))
mu <- rnorm(1,mean,sd)
alpha <- (n+nu0)/2
beta <- (Sum2+n*(mu-Xbar)^2+ramda0)/2
sig2 <- rinvgamma(1,alpha,beta)
if(i > 0) {
sample_mu[i] <- mu #μのサンプル
sample_var[i] <- sig2 #σ^2のサンプル
sample_sd[i] <- sqrt(sig2) #σのサンプル
}
}
summary(sample_mu)
summary(sample_sd)
st_mu <- c(mean(sample_mu),sqrt(var(sample_mu)))
st_sd <- c(mean(sample_sd),sqrt(var(sample_sd)))
print(st_mu)
print(st_sd)
plot(sample_mu,type="l")
par(ask=TRUE)
hist(sample_mu, breaks=seq(-0.5,0.5,0.05))
par(ask=TRUE)
plot(sample_var,type="l")
par(ask=TRUE)
hist(sample_var)

７．結果の表示
> summary(sample_mu)
Min. 1st Qu. Median Mean 3rd Qu. Max.
-0.25730 0.01666 0.05740 0.05707 0.09790 0.37190
> summary(sample_sd)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.1187 0.1719 0.1847 0.1867 0.1995 0.3342
> print(st_mu)
[1] 0.05707447 0.06054553
> print(st_sd)
[1] 0.18669853 0.02113967