千葉大学フロンティア医工学センター　羽石秀昭様　「判別分析２」.
のテキストで学習した記録です。

線形判別分析は、以下の前提条件が成立する必要があります。

各グループは多変量正規分布に従う
全てのグループが同じ分散共分散行列を持つ（等分散性）

1 判別分析の種類.

Rで判別分析いろいろ（11種類）では色々な判別方法が述べられていますが、
ここでは「線形判別分析（LDA）」「ベイズ判別」について紹介します。

1.1 フィッシャー判別分析.

HELLO CYBERNETICS 様　「フィッシャーの線形判別」（二次形式表現）によれば、
フィッシャーの線形判別は、データを分類するための都合の良い線形変換を見つける手法です。手順としては、

［手順１］群間変動最大の方向を見つける.
［手順２］誤判別確率最小となる分離点を見つける.

事前確率が等しい前提では、クラス重心の中点の射影が分離点になります。
事前確率を考慮する場合は、ベイス判別の方法を用いることになります。

［手順１］ 判別直線は、
　 \(w={S_W}^{-1}(m_1-m_2)\) を法線ベクトルとし、各クラスの平均の中点を通る直線です。

　２つの群の判別の場合、

判別直線の法線ベクトルは、２つの群が最も分離して見える方向である。
この法線方向に射影すれば、２つの群は２つの１次元正規分布になる。

　（マハラノビスから求めた判別方向は、群間変動から求めたものと一致している）

［手順２］ 次に、分離の境界点をどこにするかが問題となります。
・事前確率が同じである場合であれば、２つの正規分布関数の交点が分離の境界点となります。
（判別直線が２群の重心の中点を通る場合で、このとき誤判別確率は最小となる）

誤判別確率.
　\(P_e=T(R_A,R_B)=\pi_A\int_{R_B}f_A(x)dx+\pi_B\int_{R_A}f_B(x)dx\).
　人文2022②では、
　　　\(\displaystyle \int_{R_B}f_A(x)dx=P(Y<0\mid G_A )\)：判別関数 Y によって計算している。

・事前確率が異なる場合であれば、比率によって判別の境界が移動する。（羽石サイト2のp11）.
　ベイズ判別の方法で事後確率を比較すると、
　境界線は、羽石サイト２のp14（２群の分散共分散行列が共通のとき）

筆者は、
　　　事前確率を考慮するのならば、はじめからベイズ判別を使うべきと思う。

人文2015④(3)(4)は、2変量正規分布に対するベイズ予測ルールで、
　事前確率が等しいので重心の中点を通る判別直線になる。
　ベイズ予測ルールは事前確率が異なる場合も対応可能である。
　（問題はそこまで問うてはいないが）

詳しくは当サイト「判別分析_まとめ」を参照.

1.2 マハラノビスの距離による判別.

確率的距離が等しくなる条件を２次形式を用いて計算することになり、境界を直接算出できます。方向と分離点が同時に算出されていることになりますが、事前確率を考慮した判別はできません。

人文2022②(2)では、
マハラノビスの距離を用いて判別直線を導出しています。

1.3 ベイズ事後確率最大.

ベイズによる事後確率が最大となる群に判別する方法で、大森先生のテキストによれば、最適な判別規則ということになります。
人文2015④(3)(4)が該当する問題です。

2 分散分析と判別分析.

射影済みクラスの判別と考えることができます。判別式の係数は１になります。

library(MASS) 
x<-c(1:3);y<-c(3:5) 
d<-data.frame(x,y);d

x	y
1	3
2	4
3	5

library(reshape)
dat <- melt(d, variable_name = "group");dat

## Using  as id variables

group	value
x	1
x	2
x	3
y	3
y	4
y	5

boxplot(value ~ group,data=dat)

(Z<- lda(group~ value,data=dat))

## Call:
## lda(group ~ value, data = dat)
## 
## Prior probabilities of groups:
##   x   y 
## 0.5 0.5 
## 
## Group means:
##   value
## x     2
## y     4
## 
## Coefficients of linear discriminants:
##       LD1
## value   1

predict(Z)$class#判別の予測値

## [1] x x x x y y
## Levels: x y

round(predict(Z)$posterior,2)#判別の事後確率

##      x    y
## 1 0.98 0.02
## 2 0.88 0.12
## 3 0.50 0.50
## 4 0.50 0.50
## 5 0.12 0.88
## 6 0.02 0.98

predict(Z)$x#判別関数の値

##   LD1
## 1  -2
## 2  -1
## 3   0
## 4   0
## 5   1
## 6   2

#plot(Z,dimen=1)

判別関数は平均を通ることから、　　\(f(z)=z-3\).

3 人文2018③.

3.1 事前確率が等しいとき.

サンプル数を100ずつとする。

library(MASS) 
mA <- c(30,50)
mB <- c(40,50)
Sigma <- matrix(c(20,14,14,20), 2, 2)
Sigma

##      [,1] [,2]
## [1,]   20   14
## [2,]   14   20

dataA <- mvrnorm(100, mA, Sigma)
dataB <- mvrnorm(100, mB, Sigma)
dA <- data.frame(X=dataA[,1],Y=dataA[,2],Z=rep("A",100))
dB <- data.frame(X=dataB[,1],Y=dataB[,2],Z=rep("B",100))
plot(dA[,1:2],xlim=c(20,50),ylim=c(35,65))
par(new=T)
plot(dB[,1:2],xlim=c(20,50),ylim=c(35,65),col="red")
d <- rbind(dA,dB)
#クラス重心の中点
m1 <- c(mean(dA[,1]),mean(dA[,2]))
m2 <- c(mean(dB[,1]),mean(dB[,2]))
m <- (m1+m2)/2;m

## [1] 34.73178 49.82943

Z<- lda(Z~ .,data=d)
# ldaから求めた法線ベクトル
k <- Z $ scaling;k

##          LD1
## X  0.3016142
## Y -0.1958176

b <- m[1]*k[1]/k[2]+m[2];b

## [1] -3.667259

a <- (-1)*k[1]/k[2];a

## [1] 1.540281

abline(b,a)

3.1.1 法線ベクトルの理論値.

n <- solve(Sigma)%*%(m2-m1)#法線ベクトル
(-n[1]/n[2])

## [1] 1.391717

ldaから求めた傾きとほぼ一致している。

3.1.2 ldaの結果.

## Call:
## lda(Z ~ ., data = d)
## 
## Prior probabilities of groups:
##   A   B 
## 0.5 0.5 
## 
## Group means:
##          X        Y
## A 29.97313 50.00691
## B 39.49043 49.65196
## 
## Coefficients of linear discriminants:
##          LD1
## X  0.3016142
## Y -0.1958176

head(predict(Z)$class)#判別の予測値

## [1] A A A A A A
## Levels: A B

head(round(predict(Z)$posterior,2))#判別の事後確率

##      A    B
## 1 0.72 0.28
## 2 0.99 0.01
## 3 0.62 0.38
## 4 0.89 0.11
## 5 0.86 0.14
## 6 0.94 0.06

head(predict(Z)$x)#判別関数の値

##          LD1
## 1 -0.3281589
## 2 -1.6954015
## 3 -0.1656343
## 4 -0.7073589
## 5 -0.6304454
## 6 -0.9323848

3.2 事前確率が異なるとき.

A:B=100:1000 とする。

library(MASS) 
mA <- c(30,50)
mB <- c(40,50)
Sigma <- matrix(c(20,14,14,20), 2, 2)
Sigma

##      [,1] [,2]
## [1,]   20   14
## [2,]   14   20

dataA <- mvrnorm(1000, mA, Sigma)
dataB <- mvrnorm(100, mB, Sigma)
dA <- data.frame(X=dataA[,1],Y=dataA[,2],Z=rep("A",1000))
dB <- data.frame(X=dataB[,1],Y=dataB[,2],Z=rep("B",100))
plot(dA[,1:2],xlim=c(20,50),ylim=c(35,65))
par(new=T)
plot(dB[,1:2],xlim=c(20,50),ylim=c(35,65),col="red")
d <- rbind(dA,dB)
m1 <- c(mean(dA[,1]),mean(dA[,2]))
m2 <- c(mean(dB[,1]),mean(dB[,2]))
m <- (m1+m2)/2;m

## [1] 34.99957 49.97389

Z<- lda(Z~ .,data=d)
k <- Z $ scaling;k

##          LD1
## X  0.3229731
## Y -0.2343254

b <- m[1]*k[1]/k[2]+m[2];b

## [1] 1.733625

a <- (-1)*k[1]/k[2];a

## [1] 1.37831

abline(b,a)

n <- solve(Sigma)%*%(m2-m1)#法線ベクトル
(-n[1]/n[2])

## [1] 1.404238

LDR（線形判別）ではサンプル数の影響を受けていないようである。分離点は各クラスの重心の中点であり、サンプル数が増えても変わらないからである。

3.3 ベイズ判別.

羽石先生の資料p12-p15 を参考にして整理する。
(事後確率)＝(事前確率)×(尤度) であるから、
　　　\(P(G_A|x) \propto P(G_A)~P(x|G_A)\).
Aの判別関数　\(g_A(x)=\log P(G_A)-\dfrac{1}{2}~M_A^2\)　（\(M_A^2\) はマハラノビスの距離）.
同様に.
Bの判別関数　\(g_B(x)=\log P(G_B)-\dfrac{1}{2}~M_B^2\)　（\(M_B^2\) はマハラノビスの距離）.
２群を判別する関数は、
　　　\(f(x)=\log \dfrac{P(G_A)}{P(G_B)}-\dfrac{1}{2}(M_A^2-M_B^2)\).
　　　第２項はマハラノビスの距離の差であり、事前確率を考慮しない判別境界を表す。

\(\Sigma^{-1}=\dfrac{1}{204}\begin{pmatrix}20&-14\\-14&20\end{pmatrix}=\dfrac{1}{102}\begin{pmatrix}10&-7\\-7&10\end{pmatrix}\) であるので、

　　\(M_A^2=\dfrac{1}{102}\left[10(x_1-30)^2-14(x_1-30)(x_2-50)+10(x_2-50)^2\right]\).
　　\(M_B^2=\dfrac{1}{102}\left[10(x_1-40)^2-14(x_1-40)(x_2-50)+10(x_2-50)^2\right]\).

\(-\dfrac{1}{2}(M_A^2-M_B^2)=-\dfrac{1}{2\cdot 102}\left[200(x_1-35)-140(x_2-50)\right]=-\dfrac{70}{102}\left(\dfrac{10}{7}x_1-x_2\right)\).
これは事前確率を考慮しないときの判別直線であり、事前確率を考慮する場合は \(\log \dfrac{P(G_A)}{P(G_B)}\) だけ平行移動することになる。
　　　\(x_2=\dfrac{10}{7}x_1-\dfrac{102}{70}\log \dfrac{P(G_A)}{P(G_B)}\).

plot(dA[,1:2],xlim=c(20,50),ylim=c(35,65))
par(new=T)
plot(dB[,1:2],xlim=c(20,50),ylim=c(35,65),col="red")
abline(-102/70*log(10),10/7)

判別直線が平行移動し、Aの判別領域が広くなっている。

plot(dA[,1:2],xlim=c(20,50),ylim=c(35,65))
par(new=T)
plot(dB[,1:2],xlim=c(20,50),ylim=c(35,65),col="red")
d <- rbind(dA,dB)
Z<- lda(Z~ .,data=d)
k <- Z $ scaling;k

##          LD1
## X  0.3229731
## Y -0.2343254

m[1]*k[1]/k[2]+m[2]

## [1] 1.733625

b <- m[1]*k[1]/k[2]+m[2]-(102/70)*log((1000/100))
a <- (-1)*k[1]/k[2]
abline(b,a)

## [1] -1.621571

ldaから得られた判別直線を平行移動しても、ほぼ同じ結果が得られている。

事前確率を考慮することにより、誤判別確率を最小にすることができる。このことを計算で確かめるのが、次の人文2022②である。

4 人文2022②.

(3)は、マハラノビスの距離を利用した判別.

\(\displaystyle P(x\mid G_1 )=\dfrac{2}{2\pi\sqrt{3}}\exp\left[-\dfrac{1}{2\cdot 0.75}\left[(x_1-1)^2-(x_1-1)(x_2-1)+(x_2-1)^2\right]\right]\).

(4)は、誤判別確率最小の分離点を考える問題である。
射影によって得られる２つの正規分布関数の交点が、最適誤判別確率となる分離点である。

4.1 シミュレーション.

4.1.1 事前確率が等しいとき.

サンプル数は、100,100 とする。

library(MASS) 
mA <- c(1,1)
mB <- c(-1,-1)
Sigma <- matrix(c(1,0.5,0.5,1), 2, 2)
Sigma

##      [,1] [,2]
## [1,]  1.0  0.5
## [2,]  0.5  1.0

dataA <- mvrnorm(100, mA, Sigma)
dataB <- mvrnorm(100, mB, Sigma)
dA <- data.frame(X=dataA[,1],Y=dataA[,2],Z=rep("A",100))
dB <- data.frame(X=dataB[,1],Y=dataB[,2],Z=rep("B",100))
plot(dA[,1:2],xlim=c(-3,3),ylim=c(-3,3))
par(new=T)
plot(dB[,1:2],xlim=c(-3,3),ylim=c(-3,3),col="red")
d <- rbind(dA,dB)
mean(dA[,1]);mean(dB[,1])

## [1] 0.8781042

## [1] -0.9359348

mean(dA[,2]);mean(dB[,2])

## [1] 0.9034119

## [1] -0.9586466

#中心点
m1 <- c(mean(dA[,1]),mean(dA[,2]))
m2 <- c(mean(dB[,1]),mean(dB[,2]))
m <- (m1+m2)/2;m

## [1] -0.02891531 -0.02761734

#線形判別
Z<- lda(Z~ .,data=d)
k <- Z $ scaling;k #法線ベクトル

##          LD1
## X -0.5813098
## Y -0.5689971

b <- m[1]*k[1]/k[2]+m[2];b #切片

## [1] -0.05715835

a <- (-1)*k[1]/k[2];a #傾き

## [1] -1.021639

abline(b,a)

法線ベクトルを \(k=(a,b)\) とすると、判別直線は.
　　　\(y=-\dfrac{a}{b}x+c\).
事前確率が等しいときは平均点を通るので、
　　　\(c=\bar{x}\cdot \dfrac{a}{b}+\bar{y}\).
これが切片である。

4.1.2 事前確率が異なるとき.

サンプル数は、500,100 とする。

library(MASS) 
mA <- c(1,1)
mB <- c(-1,-1)
Sigma <- matrix(c(1,0.5,0.5,1), 2, 2)
Sigma

##      [,1] [,2]
## [1,]  1.0  0.5
## [2,]  0.5  1.0

dataA <- mvrnorm(500, mA, Sigma)
dataB <- mvrnorm(100, mB, Sigma)
dA <- data.frame(X=dataA[,1],Y=dataA[,2],Z=rep("A",500))
dB <- data.frame(X=dataB[,1],Y=dataB[,2],Z=rep("B",100))
plot(dA[,1:2],xlim=c(-3,3),ylim=c(-3,3))
par(new=T)
plot(dB[,1:2],xlim=c(-3,3),ylim=c(-3,3),col="red")
d <- rbind(dA,dB)
#中心点
m1 <- c(mean(dA[,1]),mean(dA[,2]))
m2 <- c(mean(dB[,1]),mean(dB[,2]))
m <- (m1+m2)/2;m

## [1] -0.02250538 -0.03238948

#線形判別
Z<- lda(Z~ .,data=d)
k <- Z $ scaling;k

##          LD1
## X -0.6625496
## Y -0.5076683

b <- m[1]*k[1]/k[2]+m[2];b

## [1] -0.06176088

a <- (-1)*k[1]/k[2];a

## [1] -1.305084

abline(b,a)

線形判別では判別直線に大きな変化はない。
（傾きが結構変わる印象ではある）.