「Rによるモンテカルロ法入門」読書ノート：6章 メトロポリス・ヘイスティング・アルゴリズムを用いたモデル選択

例 6.8 モデル選択

この問題の目的とデータセットの説明

線形回帰のモデル選択について考える．data(swiss)はスイスの47地域における出生率と社会経済指標のデータである．社会経済指標には，農業従事率や教育の度合い，特定の宗教の割合，幼児の死亡率など，合計で5つの標準化された観測値が含まれている．今回は，線形回帰によって出生率と社会経済指標の関係を調べることを目標としている．つまり，社会経済指標の中で何が出生率に効いているかを知りたいということだ．そのような線形回帰において一番重要なのが，5つある説明変数をどうやって組み合わせてモデルを作るかという問題である．5つの説明変数から作ることの出来るモデルは25 個なので，今回の場合のモデルの総数はたかだか32個程度だが，たとえば説明変数が20個になると，モデルの総数は220 = 1048576個となり，全部のモデルをそれぞれ計算するのは現実的ではなくなってしまう．そこで，モデルの全パターンを試すことなく，もっとも出生率に効いてる社会経済指標だけを選んでモデルを作りたい．そのために，今回は周辺分布の最大化によるモデル選択をランダム・ウォーク・メトロポリス・ヘイスティング・アルゴリズムを使って算出してみる．

数式の説明

通常の線形回帰を考えるということで，出生率yと社会経済指標xは以下の式で表される．

\bold{y} = \beta X + \varepsilon ，\varepsilon \sim \mathcal{N}_n(0,\sigma^2 I_n)

このときの\bold{y}|\beta,\sigma^2,X は \bold{y}|\beta,\sigma^2,X \sim \mathcal{N}_n(X\beta,\sigma^2I_n) となる．

ここで，出生率yは47の地域でそれぞれ観測されているため，観測数n=47のベクトルy={y_1..yn}として表される．

そして，このときの尤度関数は

l(\beta,\sigma^2|\bold{y},X)=\big(\frac{1}{\sqrt{2\pi\sigma^2}}\big)^n \exp\big[-\frac{1}{2\sigma^2}(\bold{y}-X\beta)^{\mathrm{T}}(\bold{y}-X\beta)\big]

である．

このあたりの式変形は練習問題6.5を参照．

データセットの準備

swissデータセットを読み込んだあと，応答変数yとして出生率の対数，説明変数Xとして5つの社会経済指標を代入している．説明変数には

1
2
3
4

data(swiss)

y <- log(as.vector(swiss[,1]))
X <- as.matrix(swiss[,2:6])

説明変数である社会経済指標の詳細は以下の通りである．

1
2

> names(swiss)
[1] "Fertility"        "Agriculture"      "Examination"      "Education"        "Catholic"         "Infant.Mortality"

関数の定義

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

inv <- function(X){
  EV <- eigen(X)
  EV$vector %*% diag(1/EV$values) %*% t(EV$vector)
}

lpostw <- function(gam, y, X, beta){
  n <- length(y)
  qgam <- sum(gam)
  Xt1 <- cbind(rep(1,n), X[,which(gam==1)])
  if(qgam != 0){
    P1 <- Xt1 %*% inv(t(Xt1) %*% Xt1) %*% t(Xt1)
  }else{
    P1 <- matrix(0, n, n)
  }
  -(qgam+1)/2 * log(n+1) - n/2 *
    log(t(y) %*% y
        - n/(n+1)
        * t(y) %*% P1 %*% y-1/(n+1) * t(beta)
        %*% t(cbind(rep(1,n), X)) %*% P1
        %*% cbind(rep(1,n), X) %*% beta
    )
}

gocho <- function(niter, y, X){
  lga <- dim(X)[2]
  beta <- lm(y ~ X)$coeff
  gamma <- matrix(0, nrow=niter, ncol=lga)
  gamma[1,] <- sample(c(0,1), lga, rep=T)
  for(t in 1:(niter-1)){
    j <- sample(1:lga,1)
    gam0 <- gamma[t,]
    gam1 <- gamma[t,]
    gam1[j] <- 1 - gam0[j]
    pr <- lpostw(gam0, y, X, beta)
    pr <- c(pr, lpostw(gam1, y, X, beta))
    pr <- exp(pr - max(pr))
    gamma[t+1,] <- gam0
    if(sample(c(0,1), 1, prob=pr)){
      gamma[t+1,] <- gam1
    }
  }
  gamma
}

練習問題 6.5

a. \beta|\sigma^2,X からX\beta|\sigma^2,X と\bold{y}|\sigma^2,X を導出する

1つ目は，正規分布の線形変換を使う．これは，X \sim \mathcal{N}(\mu,\sigma^2) のときY = AX + b とするとY \sim \mathcal{N}(A\mu+b,A\sigma^2 A^\mathrm{T}) となる． この性質により，

X\beta|\sigma^2,X\sim\mathcal{N}_n(X\tilde{\beta},n\sigma^2X(X^\mathrm{T}X)^{-1}X^\mathrm{T})

となる．

2つ目は，正規分布の再生性を使う．これは，X \sim \mathcal{N}(\mu_X,\sigma_X^2) とY \sim \mathcal{N}(\mu_Y,\sigma_Y^2) が独立ならばX+Y \sim \mathcal{N}(\mu_X+\mu_Y,\sigma_X^2+\sigma_Y^2) が成り立つというものである．いま，\bold{y}=\beta X+\varepsilon における\beta X と\varepsilon はともに独立の正規分布なので，1つ目で求めた式と合わせて

\bold{y}|\sigma^2,X \sim \mathcal{N}_n(X\tilde{\beta},\sigma^2(I_n+n\sigma^2X(X^\mathrm{T}X)^{-1}X^\mathrm{T}))

となる．

b. & c.

まだよく分かっていないので後回し．設問b.に関しては，どうやら「Bayesian Core: A Practical Approach to Computational Bayesian Statistics (Springer Texts in Statistics)」に詳しく書いてある(らしい)．該当部分をざっと読んでみたけど，よく分からなかった…．

参考