「Rによるモンテカルロ法入門」読書ノート：6章 メトロポリス・ヘイスティング・アルゴリズム その3

練習問題 6.13

この問題は例6.3と同じく，線形回帰におけるパラメータの事後分布からメトロポリス・ヘイスティング・アルゴリズムを用いてサンプリングを行うことを目標とする．線形回帰の対象はmcsmパッケージに収められているchallengerデータセットで，過去に打ち上げられたスペースシャトルにおけるOリングの破損の有無と，その時の気温に関するものである．

a. glm()関数を使ってロジスティック回帰を行う

打ち上げられたスペースシャトルのOリングは破損しているかしていないかの2値で表されるため，今回はロジスティック回帰を用いた当てはめを行う．Oリングが破損している場合は1を，破損していない場合には0が割り当てられており，ある気温x_i におけるOリングの破損確率は

P(x_i) = \frac{\exp(\alpha+\beta_i)}{1+\exp(\alpha+\beta_i)}

で表される．今回はこの\alpha と\beta をglm()関数を使って推定する．

なお，以下のコードではロジスティック関数を式変形して\frac{1}{1+\exp(-(\alpha+\beta_i))} として計算している．

1
2
3
4
5
6

library(mcsm)
data(challenger)

fit <- glm(oring ~ temp, data=challenger, family=binomial)
plot(challenger$temp, challenger$oring, pch=19)
curve(1/(1+exp(-(fit$coef[1]+fit$coef[2]*x))), from=53, to=81, add=T, col="orange", lwd=4)

この図は，横軸に気温（華氏），縦軸にOリングの破損の有無をプロットしたものである．また，glm()関数で推定したロジスティック曲線をオレンジ色で示しており，任意の気温における破損する確率と対応している．この図を見ると，華氏65前後で破損確率が50%となり，気温が高い場合は破損確率が低く，逆に気温が低い時には破損確率が高いことがわかる．

glm()関数のパラメータ推定の詳細は以下の通り．

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

> summary(fit)

Call:
glm(formula = oring ~ temp, family = binomial, data = challenger)

Deviance Residuals:
    Min       1Q   Median       3Q      Max
-1.0611  -0.7613  -0.3783   0.4524   2.2175

Coefficients:
            Estimate Std. Error z value Pr(>|z|)
(Intercept)  15.0429     7.3786   2.039   0.0415 *
temp         -0.2322     0.1082  -2.145   0.0320 *
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

(Dispersion parameter for binomial family taken to be 1)

    Null deviance: 28.267  on 22  degrees of freedom
Residual deviance: 20.315  on 21  degrees of freedom
AIC: 24.315

Number of Fisher Scoring iterations: 5

b. ロジスティック回帰で推定したパラメータの値を用いて，メトロポリス・ヘイスティング・アルゴリズムでサンプリングを行う

これは前回と殆ど変わらず．

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41

y <- challenger$oring
x <- challenger$temp

loglik <- function(a,b){sum(y*(a+b*x) - log(1+exp(a+b*x)))}
logq <- function(a,b){
  dnorm(a,mean=a.est[1],sd=a.std,log=T)
  + dnorm(b,mean=b.est[1],sd=b.std,log=T)
}

Niter <- 5000
n <- length(challenger$temp)
fit.param <- summary(fit)$coef
a.est <- fit.param[1,1]
a.std <- fit.param[1,2]
b.est <- fit.param[2,1]
b.std <- fit.param[2,2]


for(i in 2:Niter){
  candidate <- c(rnorm(1, mean=a.est[1], sd=a.std), rnorm(1, mean=b.est[1], sd=b.std))
  fy <- loglik(candidate[1], candidate[2])
  fx <- loglik(a.est[i-1], b.est[i-1])
  qxy <- logq(a.est[i-1], b.est[i-1])
  qyx <- logq(candidate[1], candidate[2])
  rho <- min(1, exp(fy-fx+qxy-qyx))

  if (runif(1)<rho){
    a.est[i] <- candidate[1]
    b.est[i] <- candidate[2]
  }else{
    a.est[i] <- a.est[i-1]
    b.est[i] <- b.est[i-1]
  }
}

plot(x, y, xlab="temp", ylab="oring", type="n")
for(i in (Niter-500):Niter){
  lines(x, 1/(1+exp(-(a.est[i]+b.est[i]*x))), lwd=2, col="grey")
}
curve(1/(1+exp(-(fit$coef[1]+fit$coef[2]*x))), from=53, to=81, add=T, col="orange", lwd=4)
points(x, y, pch=19)

c. 5000回のマルコフ連鎖によるサンプル列を生成し，サンプリングで得られた推定値を使ったロジスティック曲線を図示する

実際に上で書いたコードを動かして，事後分布からのサンプリングを行い，得られたサンプルのうち4500回から5000回までの反復で生成した当てはめ曲線を灰色の曲線で表している．

また，パラメータ\alpha と\beta の系列の変動，自己相関プロット，系列の分布を以下の図で示している．

1
2
3
4
5
6
7

par(mfrow=c(2,3))
plot(1:Niter, a.est, type="l", xlab="Iterations", main=expression(alpha))
acf(a.est)
plot(density(a.est))
plot(1:Niter, b.est, type="l", xlab="Iterations", main=expression(beta))
acf(b.est)
plot(density(b.est))

そして最後に，サンプリング推定値差の変動を示したのが以下の図である．glm()関数で推定した値を赤線で表示しており，そこからどれだけ離れているかを見ることができる．

1
2
3

plot(1:Niter, apply(cbind(a.est,b.est),1,function(x){loglik(x[1],x[2])+logq(x[1],x[2])}),
     type="l", xlab="Iteration", ylab="")
abline(h=(loglik(fit.param[1,1],fit.param[2,1])+logq(fit.param[1,1],fit.param[2,1])), col="red", lwd=2)

d. サンプリングで得られたパラメータの分布から，華氏60,50,40度のときの故障の確率と標準誤差を推定する

サンプリングで得られた値を用いて，ある温度における故障の確率と標準誤差を求める．なお，今回はBurn-inによる推定値の切り捨てを行なっておらず，サンプリング回数5000回で得られた推定値を全て使っている．

1
2
3
4
5
6
7

for(t in c(60,50,40)){
  p <- 1/(1+exp(-(a.est+b.est*t)))
  print(paste("temp:", t, "=>", "mean=", round(mean(p),digit=4), ", S.D.=", round(sd(p),digit=4) ))
}
[1] "temp: 60 => mean= 0.7476 , S.D.= 0.1282"
[1] "temp: 50 => mean= 0.9525 , S.D.= 0.0641"
[1] "temp: 40 => mean= 0.9895 , S.D.= 0.0248"