練習

統計

通常の回帰をベイズでやってみる。
基本のキ。

まずデータの呼び出しから。
オゾン濃度と気温の関係について（のみ）解析してみる。

data(airquality)


通常の回帰（GLMで）。

model.glm <- glm(Ozone ~ Temp, data = airquality, family = gaussian)
summary(model1.glm)

Call:
glm(formula = Ozone ~ Temp, family = gaussian, data = airquality)

Deviance Residuals: 
     Min        1Q    Median        3Q       Max  

    

  • 40.7295  -17.4086   -0.5869   11.3062  118.2705  
    • 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) -146.9955    18.2872  -8.038 9.37e-13 ***
Temp           2.4287     0.2331  10.418  < 2e-16 ***

    

  • 

      

    • -
      • 

    • 

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1 

(Dispersion parameter for gaussian family taken to be 562.3675)

    Null deviance: 125143  on 115  degrees of freedom
Residual deviance:  64110  on 114  degrees of freedom
  (37 observations deleted due to missingness)
AIC: 1067.7

Number of Fisher Scoring iterations: 2


R2WinBUGSを使ってのベイズ。

Ozone <- airquality$Ozone
Temp <- airquality$Temp
d <- list("Ozone", "Temp")

R2WinBUGSではデータをリスト化しなければならない。

para <- c("a", "b", "tau", "sigma")
inits <- c(a = 0, b = 10, tau = 100)

パラメータと初期値の設定。


モデル式をテキストファイルでRと同じディレクトリに作成。ファイル名はair.model.txt。

model
{
	for (i in 1 : 153) {
		Ozone[i] ~ dnorm(mean[i], tau)
		mean[i] <- a + b * Temp[i]
		}
	
	a ~ dnorm(0, 1.0E-5)
	b ~ dnorm(0, 1.0E-5)
	tau ~ dgamma(0.01, 0.01)
	sigma <- sqrt(1 / tau)
}

これで準備完了。


ちなみにtauが0.01のガンマ分布からとられているが，これ以上小さくすると初期値が上手く設定できない。（これが苦労した）

model.bugs <- bugs(data = d, inits = inits, parameters = para,
			model.file = "air.model.txt", debug = F,
			n.iter = 100000, n.burnin = 10000,
			bugs.directory = "C:/Program Files/WinBUGS14",
			working.directory = getwd())


とすると勝手にWinBUGSが起動。

print(model.bugs, digits = 5)

Inference for Bugs model at "air.model..txt", fit using WinBUGS,
 3 chains, each with 1e+05 iterations (first 10000 discarded), n.thin = 270
 n.sims = 1002 iterations saved
               mean       sd       2.5%        25%        50%        75%      97.5%    Rhat n.eff
a        -146.01977 19.26953 -183.76750 -159.00000 -146.55000 -132.90000 -107.42000 1.00110  1000
b           2.41756  0.24520    1.93402    2.25525    2.42700    2.58200    2.89192 1.00102  1000
tau         0.00178  0.00024    0.00133    0.00162    0.00178    0.00193    0.00231 0.99961  1000
sigma      23.85083  1.65390   20.82100   22.75250   23.72500   24.84750   27.44900 0.99961  1000
deviance 1064.83134  2.70355 1062.00000 1063.00000 1064.00000 1066.00000 1071.00000 1.00494  1000

For each parameter, n.eff is a crude measure of effective sample size,
and Rhat is the potential scale reduction factor (at convergence, Rhat=1).

DIC info (using the rule, pD = Dbar-Dhat)
pD = 3.1 and DIC = 1068.0
DIC is an estimate of expected predictive error (lower deviance is better).


model.bugs.mc <- as.mcmc(model.bugs$sims.matrix)
plot(model.bugs.mc)

で作図。