線形判別分析 LDA

• 複数のクラスがあって、それぞれのクラスの標本について、複数の変数が観察されている
• 「中心からのずれ」を定量する
  • 定量にあたって
    • 全体が１つの集団とみなし、全標本の全体の重心からのずれを定量する((SSall)
    • クラスごとに１つの集団とみなし、クラスの重心からのクラスの標本のずれを低了する(SSclass1, SSclass2,...)
    • クラスは「均一」であるとみなし、クラスの標本がすべてクラスの重心にあるものとして、それを全体とする。その上で、その全体(多数が同じところに存在している)について、全体の重心からのずれを定量する(SSbetween)
    • ただし、「ずれ」とは、という dxd行列(ただしdは変数の数)で測る
  • このとき、
  • となるから、を考えることで、以下に、クラスが全体を分けているかがわかる
• さらに、クラスが全体のずれ具合をうまく分けるにあたって、次元削減して、うまく分けている低次元だけを取り出すことも考える
• 実際、は、という行列計算とし、そうして出来た行列を固有値分解によって、「よく別ける方向」の取り出しも可能
• 何をやっているかと言うと
  • 標本全体が多次元正規分布に分布していると見て、それが楕球になっている
  • これを標準多次元正規分布に変換する線形変換が
  • この線形変換によって、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)