線形判別分析

1 イントロダクション

ロジスティック回帰分析ではY|Xの分布をモデリングした上で分類を行ったが，線形判別分析ではX|Yのモデリングから始める．それをベイズの定理でY|Xの分布に戻して判別境界を構成する．

まずはベイズの定理を簡単に復習する．ふたつの事象A, Bに対して条件付き確率P(A|B)は， \begin{align*} P(A|B) = \frac{P(A\cap B)}{P(B)} \end{align*} で定義される．ただし，P(B) > 0とする．いま，C_{1}\cap C_{2} = \phi, C_{1}\cup C_{2} = \Omegaとなる事象C_{1}, C_{2}（\Omegaは全事象）を考えると， \begin{align*} P(B) &= P(B\cap (C_{1}\cup C_{2})) = P(B\cap C_{1}) + P(B\cap C_{2}) \\ &= P(B|C_{1})P(C_1) + P(B|C_{2})P(C_{2}) \end{align*} となる．したがって， \begin{align*} P(C_{i}|B) = \frac{P(C_{i}\cap B)}{P(B)} = \frac{P(C_{i})P(B|C_{i})}{P(B|C_{1})P(C_1) + P(B|C_{2})P(C_{2})} \end{align*} を得る．これがベイズの定理である．

2 説明変数が単数の場合

2値の結果変数Y\{\in 0,1\}と説明変数X\in\mathbb{R}を持つとしよう．このとき，上記のベイズの定理において，C_{1} = \{Y=1\}, C_{2} = \{Y=0\}, B=Xと設定すると， \begin{align*} \frac{P(Y=1|X)}{P(Y=0|X)} = \frac{P(Y=1)P(X|Y=1)}{P(Y=0)P(X|Y=0)} \end{align*} となる．この値が1を超えるならばY=1と判別し，1を下回る場合はY=0と判別するのがベイズ判別ルールである．通常は対数を取って判別ルールを作る．そのため，Y=1と判別されるのは次の場合である： \begin{align*} \log\frac{P(Y=1)}{P(Y=0)} + \log P(X|Y=1) - \log P(X|Y=0) > 0 \end{align*}

ここで，具体的にX|Yをモデリングしてみる．線形判別分析では，等分散の正規分布をそれぞれ仮定する．すなわち，単変数の場合は， \begin{align*} X|Y=1 \sim N(\mu_{1},\sigma^{2}),\quad X|Y=0\sim N(\mu_{0},\sigma^{2}) \end{align*} とモデリングされる．ゆえに判別ルールは，N(\mu_{j},\sigma^{2})の密度関数をf_{j}(\cdot)とし，\pi_{j} = P(Y=j)とするとき， \begin{align*} \log\frac{\pi_{1}}{\pi_{0}} + \log f_{1}(x) - \log f_{0}(x) > 0 \end{align*} \tag{1} となる．整理して次を得る： \begin{align*} \log\frac{\pi_{1}}{\pi_{0}} + \frac{2x(\mu_{1}-\mu_{0}) - (\mu_{1}^{2} - \mu_{0}^{2})}{2\sigma^{2}} > 0 \end{align*} \tag{2} したがって，テストデータx\in \mathbb{R}が上記のルールを満たせばそれはY=1と判別される．特に\pi_{1}=\pi_{0}（1と0の割合が同じ）の場合，判別ルールは \begin{align*} x > \frac{\mu_{1} + \mu_{0}}{2} \end{align*} となり，下図の破線となる．

実際のトレーニングデータ(y_{1},x_{1}), (y_{2},x_{2}),\dots, (y_{n}, x_{n})からは次のように未知パラメータを推定する： \begin{align*} &\hat{\pi}_{j} = \frac{1}{n}\sum_{i:y_{i} = j}1,\quad j=0,1, \\ &\hat{\mu}_{j} = \frac{1}{n\hat{\pi}_{j}}\sum_{i:y_{i}=j}x_{i}, \\ &\hat{\sigma}^{2} = \frac{1}{n-2}\sum_{j=0}^{1}\sum_{i:y_{i}=j}(x_{i} - \hat{\mu}_{j})^{2}. \end{align*} つまり，これらの推定量を Equation 2 に代入したものが実際の判別ルールとなる．

では以上の単変数の場合の線形判別分析を実際に実行してみよう．まずはロジスティック回帰のときに使ったデータをもう一度用意する．

  library(ISLR)
  defa = data.frame(
    default = as.numeric(Default$default=="Yes"),
    balance = Default$balance
    )
  head(defa)

  default   balance
1       0  729.5265
2       0  817.1804
3       0 1073.5492
4       0  529.2506
5       0  785.6559
6       0  919.5885

トレーニングデータとテストデータに分ける．

  defa.tr <- defa[1:9000,]
  defa.te <- defa[9001:10000,]

まずは\hat{\pi}_{j}を計算する．

  pi1 = sum(defa.tr$default)/length(defa.tr$default)
  pi0 = 1 - pi1
  pi1;pi0

[1] 0.033

[1] 0.967

defaltする確率がとても低いデータだとわかる．次に\hat{\mu}_{j}を計算するために，Y=jに対応する説明変数をそれぞれ抽出する．

  p1 = which(defa.tr$default == 1)  #Y=1に対応するサンプル
  p0 = which(defa.tr$default == 0)  #Y=0に対応するサンプル
  x1 = defa.tr$balance[p1]; x0 = defa.tr$balance[p0]

  mu1 = mean(x1)
  mu0 = mean(x0)
  mu1;mu0

[1] 1757.044

[1] 802.2616

最後に\hat{\sigma}^{2}を計算する．

  sigma2 = (sum((x0 - mu0)^2) + sum((x1 - mu1)^2))/(dim(defa.tr)[1]-2)
  sigma2

[1] 204880.4

推定量をplug-inした判別ルールを使ってテストデータを判別してみると次の結果を得る．

  x = defa.te$balance
  log.ratio = log(pi1/pi0) + (2*x*(mu1-mu0) - (mu1^2 - mu0^2))/(2*sigma2)
  drule = as.numeric(log.ratio > 0)
  table(defa.te$default, drule)

ロジスティック回帰の判別性能とほとんど同様のようである．

3 説明変数が複数の場合

説明変数が複数，つまり\bm{X} = (X_{1}, X_{2},\dots, X_{p})^{T} \in \mathbb{R}^{p}の場合，\bm{X}|Yを多変量正規分布でモデリングする． \bm{X}|Y=1 \sim N_{p}(\bm{\mu}_{1}, \Sigma),\quad \bm{X}|Y=0 \sim N_{p}(\bm{\mu}_{0},\Sigma) ここで，\bm{\mu}_{j} = (\mu_{j1}, \mu_{j2},\dots, \mu_{jp})^{T}\in\mathbb{R}^{p}であり， \begin{align*} \Sigma = \begin{pmatrix} \sigma_{11} & \sigma_{12} & \cdots &\sigma_{1p} \\ \sigma_{21} & \sigma_{22} & \cdots &\sigma_{2p} \\ \vdots & \vdots & \ddots & \cdots \\ \sigma_{p1} & \sigma_{p2} & \cdots &\sigma_{pp} \end{pmatrix} \end{align*} はp\times pの正則な対称行列である．また，N_{p}(\bm{\mu}_{j}, \Sigma)の密度関数f_{j}(\bm{x})は\bm{x}\in\mathbb{R}^{p}に対して \begin{align*} f_{j}(\bm{x}) = \frac{1}{(2\pi)^{p/2}\sqrt{\det(\bm{\Sigma})}}\exp\bigg\{-\frac{1}{2}(\bm{x}-\bm{\mu}_{j})^{T}\bm{\Sigma}^{-1}(\bm{x}-\bm{\mu}_{j})\bigg\} \end{align*} で与えられる．これらを判別ルール Equation 1 に代入すると，テストデータ\bm{x}\in\mathbb{R}^{p}が \begin{align*} \log\frac{\pi_{1}}{\pi_{0}} +\bm{x}^{T}\bm{\Sigma}^{-1}(\bm{\mu}_{1}-\bm{\mu}_{0}) - \frac{1}{2}\bm{\mu}_{1}^{T}\bm{\Sigma}^{-1}\bm{\mu}_{1} + \frac{1}{2}\bm{\mu}_{0}^{T}\bm{\Sigma}^{-1}\bm{\mu}_{0} > 0 \end{align*} を満たすとき，Y=1と判別される．上記の判別ルールに含まれるパラメータはトレーニングデータ(y_{1}, \bm{x}_{1}),\dots, (y_{n},\bm{x}_{n})を用いて次のように推定される： \begin{gather*} \hat{\bm{\mu}}_{j}=\frac{1}{n\hat{\pi}_{j}}\sum_{i:y_{i}=j}\bm{x}_{i},\quad j=0,1 \\ \hat{\Sigma} = \frac{1}{n-2} \sum_{j=0}^{1}\sum_{i:y_{i}=j}(\bm{x}_{i} - \hat{\bm{\mu}}_{j})(\bm{x}_{i} - \hat{\bm{\mu}}_{j})^{T} \end{gather*}

これをまた実装してみよう．今度はbalanceの他にincomeの変数も用いた2次元の説明変数で判別ルールを構成する．

  defa2 = data.frame(
    default = as.numeric(Default$default=="Yes"),
    balance = Default$balance,
    income = Default$income
    )
  head(defa2)

  default   balance    income
1       0  729.5265 44361.625
2       0  817.1804 12106.135
3       0 1073.5492 31767.139
4       0  529.2506 35704.494
5       0  785.6559 38463.496
6       0  919.5885  7491.559

トレーニングデータとテストデータに分ける．

  defa2.tr <- defa2[1:9000,]
  defa2.te <- defa2[9001:10000,]

トレーニングデータにおいて散布図を描いてみる．

  p1 = which(defa2.tr$default == 1) 
  p0 = which(defa2.tr$default == 0)

  plot(x = NULL, y = NULL, xlim = c(0,max(defa2.tr$balance)+1), ylim = c(0,max(defa2.tr$income)+10), xlab = "balance", ylab = "income")
  points(defa2.tr$balance[p1], defa2.tr$income[p1], col="blue", pch=17)
  points(defa2.tr$balance[p0], defa2.tr$income[p0], pch=4)

黒がdefaultしないサンプルで，青がdefaultするサンプルである．推定量の計算は単変数の場合とほとんど同様である．

  x1 = defa2.tr[p1, -1]  #default=1となるサンプルに対応する説明変数
  x0 = defa2.tr[p0, -1]

\hat{\bm{\mu}}_{j}の計算．

  mu1 = apply(x1, 2, mean)
  mu0 = apply(x0, 2, mean)
  mu1; mu0

  balance    income 
 1757.044 31797.847

   balance     income 
  802.2616 33543.2742

\hat{\Sigma}の計算．

  dif1 = as.matrix(sweep(x1, 2, mu1, FUN = "-"))
  dif0 = as.matrix(sweep(x0, 2, mu0, FUN = "-"))
  sigma = (t(dif1) %*% dif1 + t(dif0) %*% dif0)/(dim(defa2.tr)[1] - 2)
  sigma

          balance      income
balance  204880.4   -975365.3
income  -975365.3 176809180.6

判別ルールを計算する．

  x = as.matrix(defa2.te[,-1])  #テストデータの行列
  inv.sigma = solve(sigma)  #逆行列の計算
  mu1 = as.vector(mu1)  #ベクトル化
  mu0 = as.vector(mu0)

  const = log(pi1/pi0) - c(0.5 * (t(mu1) %*% inv.sigma %*% mu1)) + c(0.5 * (t(mu0) %*% inv.sigma %*% mu0))
  log.ratio = const + x %*% (inv.sigma %*% (mu1-mu0))

  drule = as.numeric(log.ratio > 0)
  table(defa2.te$default, drule)

最後に判別境界を求める．判別境界は\bm{x}\in\mathbb{R}^{p}に対して， \log\frac{\pi_{1}}{\pi_{0}} +\bm{x}^{T}\bm{\Sigma}^{-1}(\bm{\mu}_{1}-\bm{\mu}_{0}) - \frac{1}{2}\bm{\mu}_{1}^{T}\bm{\Sigma}^{-1}\bm{\mu}_{1} + \frac{1}{2}\bm{\mu}_{0}^{T}\bm{\Sigma}^{-1}\bm{\mu}_{0} = 0 となる．上記のdefaultデータで計算してみると，\bm{x} = (x_{1}, x_{2})であるから，c = \log\pi_{1}/\pi_{0} - \bm{\mu}_{1}^{T}\bm{\Sigma}^{-1}\bm{\mu}_{1}/2 + \bm{\mu}_{0}^{T}\bm{\Sigma}^{-1}\bm{\mu}_{0}/2，\bm{\delta} = \Sigma^{-1}(\bm{\mu}_{1} - \bm{\mu}_{0}) = (\delta_{1}, \delta_{2})^{T}とおくと， c + x_{1}\delta_{1} + x_{2}\delta_{2} = 0 が判別境界となる．実際に図示すると以下のようになる．

  delta = (inv.sigma%*%(mu1-mu0))

  plot(x = NULL, y = NULL, xlim = c(0,max(defa2.tr$balance)+1), ylim = c(0,max(defa2.tr$income)+10), xlab = "balance", ylab = "income")
  points(defa2.tr$balance[p1], defa2.tr$income[p1], col="blue", pch=17)
  points(defa2.tr$balance[p0], defa2.tr$income[p0], pch=4)
  abline(a = -const/delta[2], b = -delta[1]/delta[2], col="red", lwd=3)

4 Discussion

判別結果を見てみると，defaultするサンプルを当てることが出来ていない．これは，トレーニングデータにdefaultしないサンプルが圧倒的に多く含まれているためであり，ひとつの解決策はY=0のサンプルからダウンサンプリングすることである．しかし，ここでは異なるアプローチを考えてみる．

判別ルール Equation 1 を考える．明らかに>0としている閾値部分0を小さくすると，Y=1と予測する可能性が高くなる．

  drule1 = as.numeric(log.ratio > - 1.5)
  table(defa2.te$default, drule1)

  drule2 = as.numeric(log.ratio > 1.5)
  table(defa2.te$default, drule2)

解析の目的に合わせて，この閾値をチューニングパラメータとして調整することもできる．

Exercise 5.1

5-foldクロスバリデーションを用いて，トレーニングデータのみから閾値を調整し，テストデータにおける判別性能を確認せよ．

また，閾値を色々と変更した場合，全体的な判別の性能をどう評価したらいいだろうか．その際に有用なのがROC曲線である．ROC曲線とは，パラメータを変更しながら，(FPR, TPR)をプロットしたものである．

FPRは低い方が良く，逆にTPRは高い方が良いため，ROC曲線の下側面積でその判別方法の性能評価が可能となる．

線形判別分析でROC曲線を実際に求めてみよう．

  x.tr = as.matrix(defa2.tr[,-1])  #トレーニングデータの行列
  log.ratio.tr = const + x.tr %*% (inv.sigma %*% (mu1-mu0))

  thr <- function(x){
    tmp <- table(defa2.tr$default, as.numeric(log.ratio.tr > x))
    TPR <- tmp[2,2]/(tmp[2,1] + tmp[2,2])
    FPR <- tmp[1,2]/(tmp[1,1] + tmp[1,2])
    return(c(TPR,FPR))
  }

  roc <- sapply(seq(max(log.ratio.tr)-0.0001,min(log.ratio.tr),length=1000),function(x){thr(x)})
  plot(roc[2,],roc[1,],type="l",lty=1,lwd=2,col="blue",xlim=c(0,1),ylim=c(0,1),xlab="FPR",ylab="TPR")

Exercise 5.2

ロジスティック回帰分析においても閾値をチューニングパラメータとみなしてROC曲線をトレーニングデータから描いてみよ．また線形判別分析のROC曲線と重ね書きし，どちらの性能が良いか確認せよ．