ベイズ統計 第5回演習

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

第5回R演習の内容

  • 線形回帰モデルのベイズ推定について理解を深める

パッケージMCMCpackを使った線形回帰モデルのベイズ推定

  • まずは最小二乗法により線形回帰モデルを推定してみます
    # データセット swiss を用いる
    data(swiss)
    # Fertility: 出生率、Agriculture: 農業従事男性率、Examination: 入隊試験で高得点をマークした召集兵の割合
    # Education: 小学校以上卒業した召集兵の割合、Catholic: カトリック教徒の割合、Infant.Mortality: 乳幼児死亡率
    # Fertilityを被説明変数、その他全ての変数を説明変数に用いる
    summary(lm(Fertility ~ . , data = swiss))
    # 有意な説明変数のみを採用したモデル
    swiss.lm <- lm(Fertility~Agriculture+Education+Catholic+Infant.Mortality , data=swiss)
    summary(swiss.lm)
  • パッケージMCMCpackを用いて線形回帰モデルをベイズ推定します
    library(MCMCpack)
    swiss.post <- MCMCregress(Fertility~Agriculture+Education+Catholic+Infant.Mortality, data=swiss, burnin = 1000, mcmc = 10000, thin = 1)
    summary(swiss.post)
    plot(swiss.post)
  • パッケージcodaを用いて収束判定を行います
    library(coda)
    # Gewekeの診断方法
    geweke.diag(swiss.post)
    # MCMCとburnin, thinの回数を変えてモデルを再推定します
    swiss.post2 <- MCMCregress(Fertility~Agriculture+Education+Catholic+Infant.Mortality, data=swiss, burnin = 5000, mcmc = 20000, thin = 2)
    summary(swiss.post2)
    # Gewekeの診断方法
    geweke.diag(swiss.post2)
    # Raftery-Lewisの診断方法
    raftery.diag(swiss.post2)

パッケージR2jagsを使った線形回帰モデルのベイズ推定

  • 次にRとJAGSを使ったギブス・サンプラーの計算方法を紹介します。ここではRのR2jagsパッケージを使います。
  • まず、JAGSホームページにアクセスしてJAGSをインストールします
    • Windowsの場合
      • .exeファイルをダウンロードして実行する
    • Mac OS X 10.4-10.6の場合
      • .dmgファイルをダウンロードして実行する
  • Rを起動し、R2jagsパッケージを呼び出す
    library(R2jags)
  • 次にJAGS(BUGS)ファイル(例えば、以下のサンプルファイル)をテキストファイルとして保存します(例:model1.txt)
    model{
    for(i in 1:n){
    y[i] ~ dnorm(mu[i], tau)
    mu[i] <- beta0 + beta1 * X[i,1] + beta2 * X[i,2] + beta3 * X[i,3] + beta4 * X[i,4] + beta5 * X[i,5]
    }
    beta0 ~ dnorm(0,1.0E-6)
    beta1 ~ dnorm(0,1.0E-6)
    beta2 ~ dnorm(0,1.0E-6)
    beta3 ~ dnorm(0,1.0E-6)
    beta4 ~ dnorm(0,1.0E-6)
    beta5 ~ dnorm(0,1.0E-6)
    tau ~ dgamma(0.01, 0.01)
    }
  • 作成したJAGS(BUGS)ファイルを適当な場所に置く
    • R2jagsパッケージをインストールしている場合、Rインストール先フォルダの下にある、R2jagsフォルダ下のmodelフォルダに置く
    • Windowsの場合
      • 「C:\Program Files\R\R-2.x.x\library\R2jags\model」にファイルを置く
      • R画面の「ファイル」→「ディレクトリの変更」で上記フォルダを作業領域として指定する
    • Macの場合
      • 「Macintosh HD」→「ライブラリ」→「Frameworks」→「R.frameworks」→「Resources」→「library」→「R2jags」→「model」にファイルを置く
  • R2jagsを使って計算するために、データ、初期値、パラメータ、モデルファイルを指定する
    data(swiss)
    y<-swiss[,1]
    k<-ncol(swiss)
    X<-as.matrix(swiss[,2:k])
    n<-nrow(X)
    data <- list("n", "y", "X")
    in1 <- list(beta0=0, beta1=0, beta2=0, beta3=0, beta4=0, beta5=0, tau=1)
    in2 <- list(beta0=0, beta1=1, beta2=1, beta3=1, beta4=1, beta5=1, tau=1)
    in3 <- list(beta0=0, beta1=2, beta2=2, beta3=2, beta4=2, beta5=2, tau=1)
    inits <- list(in1,in2,in3)
    parameters <- c("beta0", "beta1","beta2","beta3","beta4","beta5","tau")
    model.file <- system.file(package="R2jags", "model", "model1.txt")
  • jags()関数で計算する
    swiss.jags <- jags(data=data, inits=inits, parameters, n.iter=10000, 
    n.burnin=1000, n.chains=3, model.file=model.file)
  • 計算結果を表示する
    print(swiss.jags, digits=3)
  • 得られる結果は以下のとおり
    • mu.vect: 事後平均
    • sd.vect: 事後標準偏差
    • 95%確信区間及び四分位
    • Rhat: Gelman-Rubinの診断方法
    • DIC: 偏差情報量基準
    • pD: モデルの有効パラメータ数
  • 結果をプロットする
    plot(swiss.jags)
    traceplot(swiss.jags)
  • MCMCの収束判定
    # Gelman-Rubinの診断方法
    gelman.diag(as.mcmc(swiss.jags))
    # Gewekeの診断方法
    geweke.diag(as.mcmc(swiss.jags))
    # Raftery-Lewisの診断方法
    raftery.diag(as.mcmc(swiss.jags))
  • MCMCが収束しなかった場合、update()関数で収束させる
    swiss.fit <- update(swiss.jags)
    print(swiss.fit, digits=3)

地価データの線形回帰モデルのベイズ推定

  • 今回はこちらのデータを使います。ダウンロードして任意のフォルダに解凍してください。
  • パッケージを呼び出し、データを読み込みます
    library(R2jags) # 既にパッケージを呼び出している場合は不要
    lph <- read.table("lph.csv", sep=",", header=T)
    summary(lph)
  • 以下の要領で、JAGS用のデータ(説明変数・被説明変数・誤差項)を作成します
    # 変数の指定
    # 被説明変数
    y <- as.vector(lph$LPH)
    # 地域数
    n <- length(y)
    # 説明変数
    x <- as.matrix(cbind(
    lph$POPD,
    lph$EMP3D
    ))
    # 線形回帰モデルの誤差項(最小二乗法)
    model.lm <- lm(y~x)
    summary(model.lm)
  • 次に線形回帰モデルを推定するためのJAGSコードを作成し、lm.txtという名前でテキストファイルとして保存します
    model{
    for(i in 1:n){
    y[i] ~ dnorm(mu[i], tau)
    mu[i] <- b0+b1*x[i,1]+b2*x[i,2] 
    }
    b0 ~ dnorm(0,1.0E-6)
    b1 ~ dnorm(0,1.0E-6)
    b2 ~ dnorm(0,1.0E-6)
    tau ~ dgamma(0.01, 0.01)
    sigma <- 1/sqrt(tau)
    }
  • JAGSコードファイルが用意できたら、以下の要領でモデルを推定します
    # JAGS変数設定
    # データ
    data <- list("n", "y", "x")
    # MCMC初期値(事前情報)
    in1 <- list(b0=model.lm$coefficients[1],b1=model.lm$coefficients[2],
    b2=model.lm$coefficients[3], tau=1)
    in2 <- list(b0=model.lm$coefficients[1],b1=model.lm$coefficients[2],
    b2=model.lm$coefficients[3], tau=1)
    in3 <- list(b0=model.lm$coefficients[1],b1=model.lm$coefficients[2],
    b2=model.lm$coefficients[3], tau=1)
    inits <- list(in1,in2,in3)
    # パラメータ
    parameters <- c("b0", "b1", "b2", "tau", "sigma")
    # モデルファイル
    model.file <- system.file(package="R2jags", "model", "lm.txt")
    # MCMC
    lm.jags <- jags(data=data, inits=inits, parameters, n.iter=10000, 
    n.burnin=1000, n.chains=3, model.file=model.file)
  • 結果の出力
    print(lm.jags, digits=3)
    # jpeg(filename="fig1.jpg")
    plot(lm.jags)
    # dev.off()
    traceplot(lm.jags)
    # jpeg(filename="fig2.jpg")
    # traceplot(lm.jags, varname="b0")
    # dev.off()
    # lm.fit <- update(lm.jags)
    # print(lm.fit, digits=3)
  • MCMCの収束判定
    # Gelman-Rubinの診断方法
    gelman.diag(as.mcmc(lm.jags))
    # Gewekeの診断方法
    geweke.diag(as.mcmc(lm.jags))
    # Raftery-Lewisの診断方法
    raftery.diag(as.mcmc(lm.jags))

トップ   編集 凍結解除 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2014-07-08 (火) 14:52:11 (1955d)