ベイズ統計 第4回演習

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

第4回R演習の内容

  • マルコフ連鎖モンテカルロ法の収束判定について理解を深める
  • Geweke、Gelman-Rubin、Raftery-Lewisの方法を示す。

MCMCによる正規分布のベイズ推定

  • マルコフ連鎖モンテカルロ法(メトロポリス法)による正規母集団のサンプリング(第3回演習の内容と同じ)
    # 正規母集団の平均と分散⇒事前分布超パラメータの初期値
    mu <- 0; sig <- 1
    # 正規母集団から抽出された標本数
    n <- 30
    # 正規母集団から抽出した標本データ
    y <- rnorm(n, mu, sig)
    # 標本平均
    yvar <- mean(y)
    # 事前分布の超パラメータ
    mu0 <- 0  	# 平均
    sig0 <- 4 	# 標準偏差
    nu0 <- 1 	# 精度の事前ガンマ分布超パラメータ
    S0 <- 4 	# 同上
    # 繰り返し計算回数
    iter <- 10000
    # 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
    }

パッケージcoda

  • パッケージcodaを使う。
    library(coda)
    # MCMCへの変換
    data <- as.mcmc(posterior[1,1,])
    summary(data)
  • 事後分布の可視化:MCMCのトレース図と密度分布が表示されます。
    plot(data)
  • 時系列自己回帰を可視化する
    autocorr.plot(data)
  • M-H法の棄却率
     rejectionRate(data)

Gelman&Rubinの診断方法

  • Within-chain variance (W)とBetween-chain variance (B)を計算したのち、Potential Scale Reduction Factor (PSRF)を計算する。chain数は2つ以上あることが前提。
    thita <- rbind(posterior[1,1,], posterior[2,1,], posterior[3,1,])
    thita <- as.matrix(thita)
    n <- ncol(thita)
    k <- nrow(thita)
    B <- n*var(rowMeans(thita))
    W <- mean(apply(thita,1,"var"))
    Vhat <- (n-1)/n*W+B/n
    PSRF <- sqrt(Vhat/W)
    print(PSRF)

Gewekeの診断方法

  • Geweke統計量Z値を計算し、MCMCのburn-in期間(前半10%)と後半50%を比較し、|Z|値が十分に小さければMCMCが収束したと考えます。以下のコードはパッケージcodaのgeweke.diag()関数を加筆修正しています。
    c <- 1 # n.chain
    n <- length(posterior[c,1,])
    # MCMCの最初の10%をg1、後半50%をg2とする
    p1 <- 0.1
    p2 <- 0.5
    n1 <- ceiling(1+p1*(n-1))
    n2 <- floor(n-p2*(n-1))
    g1 <- posterior[c,1,1:(n1)]
    g2 <- posterior[c,1,n2:n]
    # g1とg2の平均
    g1.mean <- mean(g1)
    g2.mean <- mean(g2)
    # g1とg2のスペクトラム密度
    x1 <- as.matrix(g1)
    g1.spec <- numeric(ncol(x1))
    z1 <- 1:nrow(x1)
    for (i in 1:ncol(x1)) {
    	lm.out <- lm(x1[, i] ~ z1)
    	ar.out <- ar(x1[, i], aic = TRUE)
    	g1.spec[i] <- ar.out$var.pred/(1 - sum(ar.out$ar))^2
    }
    x2 <- as.matrix(g2)
    g2.spec <- numeric(ncol(x2))
    z2 <- 1:nrow(x2)
    for (i in 1:ncol(x2)) {
    	lm.out <- lm(x2[, i] ~ z2)
    	ar.out <- ar(x2[, i], aic = TRUE)
    	g2.spec[i] <- ar.out$var.pred/(1 - sum(ar.out$ar))^2
    	order[i] <- ar.out$order
    }
    # Geweke統計量Z値の計算
    g1.var <- g1.spec/length(g1)
    g2.var <- g2.spec/length(g2)
    z <- (g1.mean-g2.mean)/sqrt(g1.var+g2.var)
    print(list(z=z, p=c(p1, p2)))

Raftery&Lewisの診断方法

  • MCMCの標本期間が適切かどうかを判断します。以下のコードはパッケージcodaのraftery.diag()関数を加筆修正しています。
    # Raftery-Lewis
    data <- as.mcmc(posterior[1,1,])
    # 初期値設定
    q = 0.025 # 四分位偏差
    r = 0.005 # 推定精度
    s = 0.95 # 区間[q-r,q+r]
    converge.eps = 0.001 # 収束判定条件
    # Minimum sample size (Nmin)
    phi <- qnorm(0.5*(1+s))
    nmin <- as.integer(ceiling((q*(1-q)*phi^2)/r^2))
    # 結果出力用行列
    out <- matrix(nrow = 1, ncol = 4, 
    		dimnames = list(varnames(data, allow.null = TRUE), 
    		c("Burnin(M)", "Required sample size(N)", 
    		"Minimum sample size(Nmin)", "Dependence factor(I=(M+N)/Nmin)")))
    # 四分位偏差
    quant <- quantile(data, probs = q)
    # 四分位偏差を下回るMCMC 下回るとTRUE
    dichot <- mcmc(data <= quant, start = start(data), 
                   end = end(data), thin = thin(data))
    # Thinning intervalの設定
    kthin <- 0
    # Transision matrixの計算
    bic <- 1
    while (bic >= 0) {
    	kthin <- kthin + thin(data)
    	testres <- as.vector(window(dichot, thin = kthin))
    	testres <- factor(testres, levels = c(FALSE, TRUE))
    	newdim <- length(testres)
    	testtran <- table(testres[1:(newdim - 2)], 
    					testres[2:(newdim - 1)], testres[3:newdim])
    	testtran <- array(as.double(testtran), dim = dim(testtran))
    	g2 <- 0
    	for (i1 in 1:2) {
    		for (i2 in 1:2) {
    			for (i3 in 1:2) {
    				if (testtran[i1, i2, i3] != 0) {
    					fitted <- (sum(testtran[i1, i2, 1:2]) * 
                           sum(testtran[1:2, i2, i3]))/(sum(testtran[1:2, i2, 1:2]))
                           g2 <- g2 + testtran[i1, i2, i3] * log(testtran[i1, i2, i3]/fitted) * 2
                   }
               }
           }
      }
      bic <- g2 - log(newdim - 2) * 2
    }
    # Transision Matrix
    finaltran <- table(testres[1:(newdim - 1)], testres[2:newdim])
    alpha <- finaltran[1, 2]/(finaltran[1, 1]+finaltran[1, 2])
    beta <- finaltran[2, 1]/(finaltran[2, 1]+finaltran[2, 2])
    tempburn <- log((converge.eps * (alpha + beta))/max(alpha,beta))/(log(abs(1 - alpha - beta)))
    nburn <- as.integer(ceiling(tempburn) * kthin)
    tempprec <- ((2 - alpha - beta) * alpha * beta * phi^2)/(((alpha + beta)^3) * r^2)
    nkeep <- as.integer(ceiling(tempprec) * kthin)
    iratio <- (nburn + nkeep)/nmin
    # output matrix
    out[1, 1] <- nburn
    out[1, 2] <- nkeep + nburn
    out[1, 3] <- nmin
    out[1, 4] <- signif(iratio, digits = 3)
    # print output
    print(list(params=c(r=r, s=s, q=q), out=out))

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