ryamadaの遺伝学・遺伝統計学メモ このページをアンテナに追加 RSSフィード

数学・コンピュータ関連の姉妹ブログ『ryamadaのコンピュータ・数学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

2017-10-18 線形判別分析 LDA

[][]線形判別分析 LDA

  • 複数のクラスがあって、それぞれのクラスの標本について、複数の変数が観察されている
  • 「中心からのずれ」を定量する
    • 定量にあたって
      • 全体が1つの集団とみなし、全標本の全体の重心からのずれを定量する((SSall)
      • クラスごとに1つの集団とみなし、クラスの重心からのクラスの標本のずれを低了する(SSclass1, SSclass2,...)
      • クラスは「均一」であるとみなし、クラスの標本がすべてクラスの重心にあるものとして、それを全体とする。その上で、その全体(多数が同じところに存在している)について、全体の重心からのずれを定量する(SSbetween)
      • ただし、「ずれ」とは、(x-m) ¥cdot (x-m)^Tという dxd行列(ただしdは変数の数)で測る
    • このとき、SSall = ¥sum_{i=1}^{nclass} SSclass_i + SSbetween
    • となるから、¥frac{SSbetween}{SSall}を考えることで、以下に、クラスが全体を分けているかがわかる
  • さらに、クラスが全体のずれ具合をうまく分けるにあたって、次元削減して、うまく分けている低次元だけを取り出すことも考える
  • 実際、¥frac{SSbetween}{SSall}は、SSall^{-1}SSbetweenという行列計算とし、そうして出来た行列を固有値分解によって、「よく別ける方向」の取り出しも可能
  • 何をやっているかと言うと
    • 標本全体が多次元正規分布に分布していると見て、それが楕球になっている
    • これを標準多次元正規分布に変換する線形変換が¥SSall^{-1}
    • この線形変換によって、SSbetweenの形を変えてみると、特に伸びている方向とそうでもない方向に分かれる
    • SSbetweenが伸びている方向っていうのは、相対的にSSintraclassが小さい方向なので、この線形変換したSSbetween(SSbetween/SSall)の楕球の形を固有値分解で評価すると、SSbetweenを説明する軸が取り出せる、という仕組み
# クラスごとの標本数
n1 <- 10
n2 <- 15
n3 <- 18

# 変数の数
d <- 4
# クラスごとの標本データ
X1 <- matrix(rnorm(n1*d),ncol=d)
X2 <- matrix(rnorm(n2*d),ncol=d)
X3 <- matrix(rnorm(n3*d),ncol=d)

Xall <- rbind(X1,X2,X3)
# クラスごと、全体の変数ごとの平均
m1 <- apply(X1,2,mean)
m2 <- apply(X2,2,mean)
m3 <- apply(X3,2,mean)
mall <- apply(Xall,2,mean)

# クラスの全標本が、クラス重心にあるものとして
# 作った標本データ
X1. <- matrix(rep(m1,each=n1),ncol=d)
X2. <- matrix(rep(m2,each=n2),ncol=d)
X3. <- matrix(rep(m3,each=n3),ncol=d)

Xall. <- rbind(X1.,X2.,X3.)
ss1 <- ss2 <- ss3 <- ssbetween <- ssall <- matrix(0,d,d)
for(i in 1:n1){
	ss1 <- ss1 + matrix(X1[i,]-m1,ncol=1) %*% matrix(X1[i,]-m1,nrow=1)
}
for(i in 1:n2){
	ss2 <- ss2 + matrix(X2[i,]-m2,ncol=1) %*% matrix(X2[i,]-m2,nrow=1)
}

for(i in 1:n3){
	ss3 <- ss3 + matrix(X3[i,]-m3,ncol=1) %*% matrix(X3[i,]-m3,nrow=1)
}
for(i in 1:(n1+n2+n3)){
	ssbetween <- ssbetween + matrix(Xall.[i,]-mall,ncol=1) %*% matrix(Xall.[i,]-mall,nrow=1)
	ssall <- ssall + matrix(Xall[i,]-mall,ncol=1) %*% matrix(Xall[i,]-mall,nrow=1)
}

# 実質的にゼロ行列
ssall - ssbetween - (ss1+ss2+ss3)

ssratio <- solve(ssall) %*% ssbetween

eigen.out <- eigen(ssratio)