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

例 6.3 二次モデルの線形回帰におけるパラメータの事後分布から，メトロポリス・ヘイスティング・アルゴリズムを用いてサンプリングする

この例題では，carsのデータセットにおいて線形回帰による二次曲線のあてはめをおこない，最尤法で推定したパラメータによる事後分布から複数のサンプリングを行うことを目標にする．今回使用するデータセットは車の速度とその時の制動距離に関するデータで，Rの入門書でよく例に使われるものである．carsのデータセットに関するHelpのExampleには，本問題と似たようなプロットを描写するコード例が載っている．

この問題では，車のスピードx と制動距離y の関係について，以下式で表される二次モデルで線形回帰をおこなう．

y_{ij}=a+bx_i+cx_i^2+\varepsilon_{ij} 

1
2
3
4
5
6
7
8
9
10
11

data(cars)

y <- cars$dist
x <- cars$speed
x2 <- x^2
fit <- lm(y ~ x + x2)
s <- seq(4,25,by=0.1)
yhat <- fit$coef[1] + fit$coef[2]*s + fit$coef[3]*s^2

plot(x,y, ylim=c(0,120), xlab="speed", ylab="dist", main="data(cars)")
lines(s, yhat, lwd=4, col="red")

以下の図は，車のスピードをx軸に，制動距離をy軸に取って，carsに含まれる50個のデータをプロットし，オレンジの線で回帰直線を表示している．また，この時の線形回帰で求められた各パラメータの推定値をsummary()関数を用いて表示している．この情報から，それぞれのパラメータの推定値は

  • a = 2.47014
  • b = 0.91329
  • c = 0.09996

ということがわかる．次の練習問題 6.11では，この推定値を元に，事後分布のサンプリングをメトロポリス・ヘイスティング・アルゴリズムで行う．

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18

> summary(fit)

Call:
lm(formula = y ~ x + x2)

Residuals:
    Min      1Q  Median      3Q     Max
-28.720  -9.184  -3.188   4.628  45.152

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept)  2.47014   14.81716   0.167    0.868
x            0.91329    2.03422   0.449    0.656
x2           0.09996    0.06597   1.515    0.136

Residual standard error: 15.18 on 47 degrees of freedom
Multiple R-squared: 0.6673,  Adjusted R-squared: 0.6532
F-statistic: 47.14 on 2 and 47 DF,  p-value: 5.852e-12

練習問題 6.11 メトロポリス・ヘイスティング・アルゴリズムを用いて事後分布からのサンプリングをする（例6.3の続き）

a. P.202の図6.6を描写する

ここで\varepsilon_{ij} \sim \mathcal{N}(0,\sigma^2) とすると，この時の尤度関数は以下式のように表される．

p(y_{ij}|x_i,a,b,c,\sigma^2) = \prod_{ij} \frac{1}{\sqrt{2\pi\sigma^2}}\exp(-\frac{(y_{ij}-a-bx_i-cx_i^2)^2}{2\sigma^2})

これを整理すると．

p(y_{ij}|x_i,a,b,c,\sigma^2) \sim \bigg( \frac{1}{\sqrt{\sigma^2}} \bigg)^N \exp{(-\frac{-1}{2\sigma^2} \sum_{ij}(y_{ij}-a-bx_i-cx_{ij})^2 )}

となる．この尤度関数をx_i,a,b,c,\sigma^2 による事後分布とみなす．それぞれのパラメータの分布は，先ほどの線形回帰で求めた推定値と標準誤差を用いる．

  • a \sim \mathcal{N}(2.47, (14.8)^2)
  • b \sim \mathcal{N}(0.913, (2.03)^2)
  • c \sim \mathcal{N}(0.100, (1.52)^2)
  • \sigma^{-2} \sim \mathcal{G}(n/2, (n-3)(15.18)^2)

なお，上の推定値や標準誤差は本書に例示されている値と若干異なっているが，今回は自分の計算環境で計算した推定値を用いることにした．

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
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60

# 例6.3での線形回帰
data(cars)
y <- cars$dist
x <- cars$speed
x2 <- x^2
fit <- lm(y ~ x + x2)

loglik <- function(a,b,c,sigma){-(n/2)*log(sigma)-sum((y-a-b*x-c*x2)^2)/(2*sigma)}
logq <- function(a,b,c,sigma){
  dnorm(a,mean=a.est[1],sd=a.std,log=T)
  + dnorm(b,mean=b.est[1],sd=b.std,log=T)
  + dnorm(c,mean=c.est[1],sd=c.std,log=T)
  - (n/2)*log(sigma) - (n-3)*(15.18)^2/(2*sigma)
}

Niter <- 4000
n <- length(x)
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]
c.est <- fit.param[3,1]
c.std <- fit.param[3,2]
sig <- 1/rgamma(1, n/2, rate=(n-3)*(15.18)^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),
                 rnorm(1, mean=c.est[1], sd=c.std),
                 1/rgamma(1, n/2, rate=(n-3)*(15.18)^2/2)
                 )
  fy <- loglik(candidate[1], candidate[2], candidate[3], candidate[4])
  fx <- loglik(a.est[i-1], b.est[i-1], c.est[i-1], sig[i-1])
  qxy <- logq(a.est[i-1], b.est[i-1], c.est[i-1], sig[i-1])
  qyx <- logq(candidate[1], candidate[2], candidate[3], candidate[4])
  rho <- min(1, exp(fy-fx+qxy-qyx))

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

s <- seq(4,25,by=0.1)
yhat <- fit$coef[1] + fit$coef[2]*s + fit$coef[3]*s^2

plot(x,y, ylim=c(0,120), xlab="speed", ylab="dist", main="data(cars)", type="n")
for(i in (Niter-500):Niter){
  lines(x, a.est[i]+b.est[i]*x+c.est[i]*x2, lwd=2, col="grey")
}
lines(s, yhat, lwd=4, col="orange")
points(x, y, pch=19)

以下の図は，先ほどの線形回帰のプロットに加えて，サンプリングから得られた値を用いて回帰曲線を灰色の曲線として描写したものである．

b. 各パラメータの収束をモニタリングする

先ほどのサンプリング結果を，a,b,c,\sigma^2 について，4000回のシミュレーションで生成された系列と，そのときの自己相関プロット，そして系列の分布をプロットした．

1
2
3
4
5
6
7
8
9
10
11
12
13

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

c. 95%信用区間を求める

1
2
3

> quantile(a.est, c(0.025,0.975))
     2.5%     97.5%
-14.90933  18.98784