ESTRELA連載 第81回

今回使うデータ

JAGSの準備

  • R2jagsを参照してください。

データの作成と表示

  • パッケージの読み込み
    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)
    }
  • Rコード
    # 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)

トップ   編集 凍結解除 差分 バックアップ 添付 複製 名前変更 リロード   新規 一覧 単語検索 最終更新   ヘルプ   最終更新のRSS
Last-modified: 2010-09-21 (火) 22:56:57 (3226d)