演習の概要

  • マルコフ連鎖モンテカルロ法による事後分布のベイズ推定について理解を深める。
  • Metropolis-Hasting法、Gibbs-sampler、正規分布のギブス・サンプリングについて演習する。

注意事項

  • 慶應義塾大学SFCで開講している「ベイズ統計」の授業履修者向けの演習用ページです。
  • 必ずしも全てのバージョンのRやOSで動作確認を行っていません。
  • Rの更新などにより、Rコードなどを予告無しに変更する場合があります。

最尤法とベイズ法

正規分布の対数尤度関数と最尤法

# 正規乱数を生成させる個数
n <- 50
# 正規乱数を生成
y <- rnorm(n, 0, 1)
# 対数尤度関数
ll <- function(x){
    mu <- x[1]; sig <- sqrt(x[2])
    sum(dnorm(y, mu, sig, log=T))   # log=T で対数化
}
# 対数尤度関数を最大化する。 
ll.max <- optim(c(0.1,0.9),ll,control=list(fnscale=-1, ndeps=1e-5))
#尤度が最大となる点
ll.max$par
## [1] 0.07473216 0.87319606
xvar <- seq(-1,1,length=50)
sig2 <- seq(0.5,2,length=50)
# 尤度の計算
llz <- matrix(0,nrow=length(xvar),ncol=length(sig2))
 for(i in 1:length(xvar)){
    for(j in 1:length(sig2)){
        llz[i,j]<-f(c(xvar[i],sig2[j]))
    }
}
## Error in f(c(xvar[i], sig2[j])): could not find function "f"
# 2次元上に結果を表示
# x座標=平均, y座標=分散
contour(xvar,sig2,llz)
# 最尤解
points(ll.max$par[1],  ll.max$par[2],pch="+",col="red")

  • 正規分布の平均を最尤法により点推定する
mu <- 0; sigma <- 1
n <- 50 # データ数
x <- rnorm(n, mu, sigma)
# 対数尤度関数
ll <- function(x, muhat, sigmahat){
 sum(dnorm(x, muhat, sigmahat, log=T))
}
# 尤度関数の値
llval <- function(x){return(ll(x, mu, sigma))}
llf <- lapply(seq(-5, 5, by=0.5), llval)
plot (seq(-5, 5, by=0.5), llf , type="o", xlab = "Values of parameter", ylab = "Log likelihood", lwd=2)

MCMCによる正規分布の最尤推定

  • マルコフ連鎖モンテカルロ法(メトロポリス法)による正規母集団のサンプリング
# 正規母集団の平均と分散⇒事前分布超パラメータの初期値
mu <- 0; sig <- 1
# 正規母集団から抽出された標本数
n <- 30
# 正規母集団から抽出した標本データ
y <- rnorm(n, mu, sig)
# 標本平均
yvar <- mean(y)
# 事前分布の超パラメータ
mu0 <- 0    # 平均
sig0 <- 4   # 標準偏差
nu0 <- 1    # 精度の事前ガンマ分布超パラメータ
S0 <- 4     # 同上
# 繰り返し計算回数
iter <- 5000
# MCMCのチェーン数
n.chains <- 3
# 事後分布
posterior <- array(dim=c(n.chains,2,iter))
# 受容される分布
accepted <- array(dim=c(n.chains, iter-1))
# メトロポリス法
for(c in 1:n.chains){
    # 事後分布
    theta.post <- array(dim=c(2, iter))
    # 初期値
    mu.init <- rnorm(1,mu0,sqrt(sig0))
    #tau.init <- rgamma(1,nu0/2,S0/2)
    #sig.init <- sqrt(1/tau.init)
    sig.init <- runif(1, 0, 2)
    theta.post[1,1] <- mu.init
    theta.post[2,1] <- sig.init
    for(t in 2:iter){
        # 需要棄却法
        # 提案分布
        h.theta <- c(theta.post[1,t-1]+rnorm(1,0,0.1),
                     theta.post[2,t-1]+rnorm(1,0,0.1))
        # 事後分布(対数)=対数尤度+事前平均(対数)+事前分散(対数)
        p.post <- sum(dnorm(y, h.theta[1], h.theta[2], log=T)) 
                    + dnorm(h.theta[1], mu0, sig0, log=T) 
                    + dunif(h.theta[2], 0, 2, log=T)
        p.prev <- sum(dnorm(y, theta.post[1,t-1], theta.post[2,t-1], log=T)) 
                    + dnorm(theta.post[1,t-1], mu0, sig0, log=T) 
                    + dunif(theta.post[2,t-1], 0, 2, log=T)
        lr <- p.post - p.prev
        r <- exp(lr)
        accept <- rbinom(1,1,prob=min(r,1))
        accepted[c,t-1] <- accept
        if(accept==1){
            theta.post[,t] <- h.theta
        }
        else {
            theta.post[,t] <- c(theta.post[1,t-1], theta.post[2,t-1])
        }
    }
    posterior[c,,] <- theta.post
}
# 事後分布(5000回)の可視化
plot(seq(1:5000),posterior[1, 1, 1:5000], type="l", xlab="Draws", ylab="Posterior", ylim=c(-4,4))

正規分布のギブス・サンプリング

  • 二変量正規分布のギブス・サンプリング
T<-n<-1000; ro<-0.5; y0<--4; mu1<-mu2<-0; sig1<-sig2<-1 
y1 <- rep(y0,n); y2<-rep(y0,n) 
for(i in 2:T){ 
y2[i] <- rnorm(1, (mu2+ro*sig2/sig1*(y1[i-1]-mu1)), sqrt(sig2*(1-ro^2))) 
y1[i] <- rnorm(1, (mu1+ro*sig1/sig2*(y2[i-1]-mu1)), sqrt(sig1*(1-ro^2)))} 
plot(y1,y2,type="o",xlim=c(-10, 10),ylim=c(-10, 10)) 

require(sm)
par(mfrow=c(1, 2))
sm.density(x=cbind(y1[1:T], y2[1:T]),
           xlab="μ1", ylab="μ2",
           zlab="", zlim=c(0, 0.15),xlim=c(-4, 4), ylim=c(-4, 4), col="white")
## Warning: weights overwritten by binning
## Loading required package: rpanel
## Warning in library(package, lib.loc = lib.loc, character.only = TRUE,
## logical.return = TRUE, : there is no package called 'rpanel'
plot(y1,y2,type="o",xlim=c(-5, 5),ylim=c(-5, 5))

# 2つの正規分布
par(mfrow=c(2, 1))
plot(seq(-4,4,0.1), dnorm(seq(-4,4,0.1),mean=mean(y1), sd=sd(y1)), xlab="", ylab="",
type="l", main="y1")
plot(seq(-4,4,0.1), dnorm(seq(-4,4,0.1),mean=mean(y2), sd=sd(y2)), xlab="", ylab="",
type="l", main="y2")

# サンプリングの結果
par(mfrow=c(1, 1))
plot(seq(1:1000), y1[1:T], type="l", xlab="Draws", ylab="Posterior", ylim=c(-4,4))
abline(-0.6,c(0,1000), col="red", lty=2, lwd=2)
abline(0.6,c(0,1000), col="red", lty=2, lwd=2)

  • 正規母集団のギブス・サンプリング
# 正規母集団の平均と分散⇒事前分布超パラメータの初期値
mu <- 0; sig <- 1
# 正規母集団から抽出された標本数
n <- 10
# 正規母集団から抽出した標本データ
y <- rnorm(n, mu, sig)
# 標本平均
yvar <- mean(y)
# 事前分布の超パラメータ
mu0 <- 0   # 平均
sig0 <- 4 # 標準偏差
nu0 <- 1 # 精度の事前ガンマ分布超パラメータ
S0 <- 4  # 同上
# 繰り返し計算回数
iter <- 1000
# MCMCのチェーン数
n.chains <- 3
# 事後分布
posterior <- array(dim=c(n.chains, 2, iter))
# ギブス・サンプリング
for(c in 1:n.chains){
    # 事後分布
    theta.post <- array(dim=c(2, iter))
    # 初期値
    mu.init <- rnorm(1,mu0,sqrt(sig0))
    tau.init <- rgamma(1, nu0/2, S0/2)
    sig.init <- sqrt(1/tau.init)
    nu.init <- runif(1, 0, 1)
    S1.init <- runif(1, 0, 1)
    theta.post[1,1] <- mu.init
    theta.post[2,1] <- sig.init
    for(t in 2:iter){
               # 平均
        mu1 <- (((mu0/sig0^2)+(n*yvar/theta.post[2,t-1]^2))
                /((1/sig0^2)+(n/theta.post[2,t-1]^2)))
               # 精度
        tau1 <- (1/sig0^2)+(n/theta.post[2,t-1]^2)
               # 標準偏差
        sig1 <- sqrt(1/tau1)
               # (逆)ガンマ分布のスケールパラメータ
        nu1 <- nu.init+n
        S1 <- S0+sum((y-theta.post[1,t-1])^2) 
        # S1 <- (nu0*S0+sum((y[i]-mu.init)^2))/nu1
               # 事後情報
        mu.post <- rnorm(1,mu1,sig1)
        tau.post <- rgamma(1,nu1/2,S1/2)
        sig.post <- sqrt(1/tau.post)
        theta.post[,t]<-c(mu.post, sig.post)
    }
    posterior[c,,] <- theta.post
}
# MCMCの結果の表示
plot(seq(1:1000),posterior[1, 1, 1:1000], type="l", xlab="Draws", ylab="Posterior",   ylim=c(-20,20), col="black")
lines(seq(1:1000),posterior[2, 1, 1:1000], col="blue")
lines(seq(1:1000),posterior[3, 1, 1:1000], col="green")
# 初期値
points(1,posterior[1, 1, 1], cex=1.5, pch=16, col="black")
points(1,posterior[2, 1, 1], cex=1.5, pch=16, col="blue")
points(1,posterior[3, 1, 1], cex=1.5, pch=16, col="green")