ベイズ統計R演習3

注意事項

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

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

• 正規分布の最尤解を求める方法 (参考: http://lbm.ab.a.u-tokyo.ac.jp/~omori/meiji2/sec5/sec5.html)

# 正規乱数を生成させる個数
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))

# ２つの正規分布
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")