ベイズ統計 第3回演習

  • ここでは、慶應大学SFCでの開講科目「ベイズ統計」での演習で用いている分析例やRコード・JAGSコードなどを、授業履修者のために示しています。
  • 必ずしも全てのバージョンのRやOSで動作確認を行っていません。
  • Rパッケージやデータの更新などにより、本書執筆時とはデータやRコードを予告無しに変更する場合があります。 本書記載のR分析例以外は、R最新版等で動作確認できていないRコードも存在しますが、参考のために記載しておきます(旧バージョンで動作することがあります)。

第3回R演習の内容

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

最尤法とベイズ法の違い

  • 最尤法は尤度関数を最大化する事により未知パラメータを点推定する方法。例えば、正規分布の最尤法は以下のようになる。(参考URL>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
    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]))
    	}
    }
    # 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による正規分布のベイズ推定

  • マルコフ連鎖モンテカルロ法(メトロポリス法)による正規母集団のサンプリング:参考URL
    # 正規母集団の平均と分散⇒事前分布超パラメータの初期値
    mu <- 0; sig <- 1
    # 正規母集団から抽出された標本数
    n <- 30
    # 正規母集団から抽出した標本データ
    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))
    # 受容される分布
    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
    }
  • 事後分布の可視化
    # 事後分布(1000回)の可視化
    plot(seq(1:1000),posterior[1, 1, 1:1000], type="l", xlab="Draws", ylab="Posterior", ylim=c(-4,4))
    # 事後分布(5000回)の可視化
    plot(seq(1:5000),posterior[1, 1, 1:5000], type="l", xlab="Draws", ylab="Posterior", ylim=c(-4,4))
    # 事後分布(10000回)の可視化
    plot(seq(1:10000),posterior[1, 1, 1:10000]
    , 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")
    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")
    # サンプリングの結果
    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")

トップ   編集 凍結解除 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2014-06-23 (月) 10:05:36 (1944d)