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

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

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

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] 35.18948 50.22268

判別関数の手計算.

#判別関数の手計算
d[1:6,1:2]
X Y
31.71100 51.90537
39.95168 61.00549
38.00671 57.83068
29.93839 53.62471
29.50599 47.61580
25.51773 50.61455 
#群内分散
var(d[1:100,1:2])
##          X        Y
## X 22.74451 15.81835
## Y 15.81835 19.26158
var(d[101:200,1:2])
##          X        Y
## X 18.71895 14.33704
## Y 14.33704 22.93239
#プールした群内分散
pSigma <- (var(d[1:100,1:2])+var(d[101:200,1:2]))*99/198
p1 <- mahalanobis(d[,1:2],m1,pSigma)
P1 <-exp((-1/2)*p1)/ (2*pi*sqrt(det(pSigma)));#P1
p2 <- mahalanobis(d[,1:2],m2,pSigma)
P2 <- exp((-1/2)*p2)/(2*pi*sqrt(det(pSigma)));#P2
P <- log(P1)-log(P2)

P <- log(P1/P2)-mean(log(P1/P2))
mean(P)
## [1] 4.440892e-17
mean(log(P1/P2))
## [1] 1.083578e-15
var(P[1:100])
## [1] 9.737009
var(P[101:200])
## [1] 9.620394
#判別関数の群内分散が１になるようにするとIDRの結果と一致する
varp <- (var(P[1:100])+var(P[101:200]))*99/198
varp
## [1] 9.678701
par(mfrow=c(3,1))
hist(P[1:100]/sqrt(varp),xlim=c(-15,15))
hist(P[101:200]/sqrt(varp),xlim=c(-15,15))
hist(P[1:200]/sqrt(varp),xlim=c(-15,15))

var(P/sqrt(varp))
## [1] 3.426809
head(P/sqrt(varp))
## [1] 1.4777231 0.8966895 0.8046911 2.4230519 1.2194094 3.1523059

判別関数の群内分散が１になるようにするとIDRの結果と一致する.

LDAによる計算.

(Z<- lda(Z~ .,data=d))
## Call:
## lda(Z ~ ., data = d)
## 
## Prior probabilities of groups:
##   A   B 
## 0.5 0.5 
## 
## Group means:
##          X        Y
## A 30.20301 50.11230
## B 40.17595 50.33306
## 
## Coefficients of linear discriminants:
##          LD1
## X  0.3168891
## Y -0.2231124
# ldaから求めた法線ベクトル
k <- Z $ scaling;k
##          LD1
## X  0.3168891
## Y -0.2231124
b <- m[1]*k[1]/k[2]+m[2];b
## [1] 0.2426448
a <- (-1)*k[1]/k[2];a
## [1] 1.420312
plot(Z,dimen=1)
abline(b,a)

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

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

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

3.1.2 ldaの結果.

Z
## Call:
## lda(Z ~ ., data = d)
## 
## Prior probabilities of groups:
##   A   B 
## 0.5 0.5 
## 
## Group means:
##          X        Y
## A 30.20301 50.11230
## B 40.17595 50.33306
## 
## Coefficients of linear discriminants:
##          LD1
## X  0.3168891
## Y -0.2231124
head(predict(Z)$class)#判別の予測値
## [1] A A A A A A
## Levels: A B
head(round(predict(Z)$posterior,2))#判別の事後確率
##      A    B
## 1 0.99 0.01
## 2 0.94 0.06
## 3 0.92 0.08
## 4 1.00 0.00
## 5 0.98 0.02
## 6 1.00 0.00
#判別関数の値
head(predict(Z)$x)
##          LD1
## 1 -1.4777231
## 2 -0.8966895
## 3 -0.8046911
## 4 -2.4230519
## 5 -1.2194094
## 6 -3.1523059
var(predict(Z)$x)
##          LD1
## LD1 3.426809
library("dplyr")
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:reshape':
## 
##     rename
## The following object is masked from 'package:MASS':
## 
##     select
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
predict <- as.data.frame(predict(Z)$x)
head(predict)
LD1
-1.4777231
-0.8966895
-0.8046911
-2.4230519
-1.2194094
-3.1523059 
dat <- mutate(predict,Z=c(rep(-1,100),rep(1,100)))
head(dat)
LD1 Z
-1.4777231 -1
-0.8966895 -1
-0.8046911 -1
-2.4230519 -1
-1.2194094 -1
-3.1523059 -1 
Z<- lda(Z~ .,data=dat)
plot(Z,dimen=2)
abline(0,1)

var(predict[1:100,])
## [1] 1.006024
var(predict[101:200,])
## [1] 0.9939757
#プール分散
(var(predict[1:100,])*99+var(predict[101:200,])*99)/200
## [1] 0.99
mA <- mean(predict[1:100,]);mA
## [1] -1.555531
mB <- mean(predict[101:200,]);mB
## [1] 1.555531
predict[1:6,]
## [1] -1.4777231 -0.8966895 -0.8046911 -2.4230519 -1.2194094 -3.1523059
P1 <- dnorm(predict[1:6,],mA,1);P1
## [1] 0.3977365 0.3211091 0.3009478 0.2738335 0.3770312 0.1114940
P2 <- dnorm(predict[1:6,],mB,1);P2
## [1] 4.008843e-03 1.972970e-02 2.461839e-02 1.457675e-04 8.488147e-03
## [6] 6.139295e-06
#事後確率
P1/(P1+P2)
## [1] 0.9900214 0.9421143 0.9243829 0.9994680 0.9779826 0.9999449

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

A:B=1000: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(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.94349 50.15914
(Z<- lda(Z~ .,data=d))
## Call:
## lda(Z ~ ., data = d)
## 
## Prior probabilities of groups:
##          A          B 
## 0.90909091 0.09090909 
## 
## Group means:
##          X        Y
## A 29.86527 49.86680
## B 40.02171 50.45149
## 
## Coefficients of linear discriminants:
##          LD1
## X  0.3104140
## Y -0.2012189
k <- Z $ scaling;k
##          LD1
## X  0.3104140
## Y -0.2012189
b <- m[1]*k[1]/k[2]+m[2];b
## [1] -3.747055
a <- (-1)*k[1]/k[2];a
## [1] 1.542668
abline(b,a)

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

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

判別関数の手計算.

var(d[1:1000,1:2])
##          X        Y
## X 19.61238 13.55002
## Y 13.55002 19.71111
var(d[1001:1100,1:2])
##          X        Y
## X 17.99292 11.88159
## Y 11.88159 19.74013
#プールした群内分散
pSigma <- (var(d[1:1000,1:2])*999+var(d[101:200,1:2])*99)/1098
p1 <- mahalanobis(d[,1:2],m1,pSigma)
P1 <-exp((-1/2)*p1)/ (2*pi*sqrt(det(pSigma)));#P1
p2 <- mahalanobis(d[,1:2],m2,pSigma)
P2 <- exp((-1/2)*p2)/(2*pi*sqrt(det(pSigma)));#P2

判別関数.

#判別関数
mahalanobis(m,m1,pSigma)
## [1] 2.317319
mahalanobis(m,m2,pSigma)
## [1] 2.317319
#事前確率による修正
P <- log(P1/P2)-mean(log(P1/P2))
mean(P)
## [1] -1.589032e-15
mean(log(P1/P2))
## [1] 3.791977
log(1000/100)
## [1] 2.302585
var(P[1:1000])
## [1] 9.281404
var(P[1001:1100])
## [1] 9.795843
#判別関数の群内分散が１になるようにするとIDRの結果と一致する
varp <- (var(P[1:1000])*999+var(P[1001:1100])*99)/1098
varp
## [1] 9.327788
par(mfrow=c(3,1))
hist(P[1:1000]/sqrt(varp),xlim=c(-15,15))
hist(P[1001:1100]/sqrt(varp),xlim=c(-15,15))
hist(P[1:1100]/sqrt(varp),xlim=c(-15,15))

head(P/sqrt(varp))
## [1]  1.2854222 -0.9265079  0.2679947 -0.2505863  0.9204135  0.6037374
var(P[1:1000]/sqrt(varp))
## [1] 0.9950273
var(P[1001:1100]/sqrt(varp))
## [1] 1.050179
#判別関数の値
head(predict(Z)$x)
##          LD1
## 1 -1.2858711
## 2  0.9250043
## 3 -0.2719654
## 4  0.2462930
## 5 -0.9169894
## 6 -0.5919564
mean(predict(Z)$x)
## [1] 8.114721e-16
mean(P/sqrt(varp))
## [1] 3.22974e-17
par(mfrow=c(3,1))
hist(log(1000/100)+log(P1/P2))
hist(log(P1/P2))

plot(Z,dimen=2)

par(mfrow=c(2,1))
hist(P[1001:1100]/sqrt(varp),xlim=c(-5,5))
hist(P[1:1000]/sqrt(varp),xlim=c(-5,5))

head(round(predict(Z)$posterior,5))
##         A       B
## 1 0.99995 0.00005
## 2 0.96316 0.03684
## 3 0.99899 0.00101
## 4 0.99515 0.00485
## 5 0.99986 0.00014
## 6 0.99962 0.00038
head(1000*P1/(1000*P1+100*P2))
## [1] 0.9999555 0.9632040 0.9990063 0.9951756 0.9998644 0.9996434

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))
## Call:
## lda(Z ~ ., data = d)
## 
## Prior probabilities of groups:
##          A          B 
## 0.90909091 0.09090909 
## 
## Group means:
##          X        Y
## A 29.86527 49.86680
## B 40.02171 50.45149
## 
## Coefficients of linear discriminants:
##          LD1
## X  0.3104140
## Y -0.2012189
k <- Z $ scaling;k
##          LD1
## X  0.3104140
## Y -0.2012189
m[1]*k[1]/k[2]+m[2]
## [1] -3.747055
b <- m[1]*k[1]/k[2]+m[2]-(102/70)*log((1000/100))
a <- (-1)*k[1]/k[2]
abline(b,a)

b
## [1] -7.10225

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

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

4 人文2022②.

判別関数.
　　　\(P(A|x)>P(B|x) \Leftrightarrow \pi_Af_A(x)>\pi_Bf_B(x)\Leftrightarrow \dfrac{\pi_Af_A(x)}{\pi_Bf_B(x)}>1\)と考え、.

判別関数\(f_q(x)\)として、
　　　\(f_q(x)=\log\left[\dfrac{\pi_Af_A(x)}{\pi_Bf_B(x)}\right]=\log\pi_A+\log f_A(x)-\log\pi_B-\log f_B(x)\).
とするなら、

\(X\) が部分母集団 Aに由来と判定する条件は、\(f_q(x)>0\)　と表される。.

特に、
部分母集団が正規分布で記述されているときは，簡単な判別規則をつくることができる．
平均 μ，分散共分散行列 Σ の ２次元正規分布の密度関数は，
　　　\(f(x,y)=\dfrac{1}{\left(\sqrt{2\pi}\right)^2\sqrt{|\Sigma|}}\exp\left[{-\dfrac{1}{2}D^2}\right] = \dfrac{1}{\left(\sqrt{2\pi}\right)^2\sqrt{|\Sigma|}}\exp\left[{-\dfrac{1}{2}D^2}\right]\).

\(P(y=1|x)~,~P(y=-1|x)\) の大小を比較するため、\(\left. P(y=1|x) \right/ P(y=-1|x)\) と1を比べる。

これは、\(f_q(x)=\log\left[\left. P(y=1|x) \right/ P(y=-1|x)\right]\) と０を比べることと同じである。

　\(P(y=1|x)=\dfrac{P(y=1)P(x|y=1)}{P(x)}\)、\(P(y=-1|x)=\dfrac{P(y=-1)P(x|y=-1)}{P(x)}\)　より、
　　 　\(\dfrac{P(y=1|x)}{P(y=-1|x)}=\dfrac{P(y=1)P(x|y=1)}{P(y=-1)P(x|y=-1)}=\dfrac{\pi_1P(x|y=1)}{\pi_2P(x|y=-1)}\).

よって、
　　　　\(f_q(x)=\log{\dfrac{\pi_1}{\pi_2}}+\log{\dfrac{P(x|y=1)}{P(x|y=-1)}}\)　　　

  • \(P(x|y=1)=\dfrac{1}{\left(\sqrt{2\pi}\right)^2\sqrt{|\Sigma_1|}}\exp\left[{-\dfrac{1}{2}(x-\mu_1)'\ \Sigma_1^{-1}(x-\mu_1)}\right]\)
  • \(P(x|y=-1)=\dfrac{1}{\left(\sqrt{2\pi}\right)^2\sqrt{|\Sigma_2|}}\exp\left[{-\dfrac{1}{2}(x-\mu_2)'\ \Sigma_2^{-1}(x-\mu_2)}\right]\) 　より、

　　\(f_q(x)=\log{\dfrac{\pi_1}{\pi_2}}+\log{\dfrac{\sqrt{|\Sigma_2|}}{\sqrt{|\Sigma_1|}}}-\dfrac{1}{2}(x-\mu_1)'\ \Sigma_1^{-1}(x-\mu_1)+\dfrac{1}{2}(x-\mu_2)'\ \Sigma_2^{-1}(x-\mu_2)\).

正規母集団で事前確率が等しく、かつ分散共分散行列が等しければ、

　　　\(f_q(x)=-\dfrac{1}{2}D_1^2+\dfrac{1}{2}D_2^2\).

これが判別関数である。

\(D_1^2=\begin{pmatrix}x_1-1&x_2-1\end{pmatrix}\Sigma^{-1}\begin{pmatrix}x_1-1\\x_2-1\end{pmatrix}=\dfrac{1}{0.75}\left[(x_1-1)^2-(x_1-1)(x_2-1)+(x_2-1)^2\right]\).

\(D_2^2=\begin{pmatrix}x_1+1&x_2+1\end{pmatrix}\Sigma^{-1}\begin{pmatrix}x_1+1\\x_2+1\end{pmatrix}=\dfrac{1}{0.75}\left[(x_1+1)^2-(x_1+1)(x_2+1)+(x_2+1)^2\right]\).
辺々引くと、

　　　\(D_1^2-D_2^2=\dfrac{1}{0.75}\left[-4x_1-4x_2+2(x_1+x_2)\right]=\dfrac{1}{0.75}\left[-2x_1-2x_2\right]\).
判別関数は. 　　　\(f_q(x)=-\dfrac{1}{2}D_1^2+\dfrac{1}{2}D_2^2=\dfrac{1}{0.75}\left[x_1+x_2\right]\).

本問では (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)は、誤判別確率最小の分離点を考える問題である。
射影によって得られる２つの正規分布関数の交点が、最適誤判別確率となる分離点である。
　　\(f_q(x)=\log{\dfrac{\pi_1}{\pi_2}}+\log{\dfrac{\sqrt{|\Sigma_2|}}{\sqrt{|\Sigma_1|}}}-\dfrac{1}{2}(x-\mu_1)'\ \Sigma_1^{-1}(x-\mu_1)+\dfrac{1}{2}(x-\mu_2)'\ \Sigma_2^{-1}(x-\mu_2)\).
において、分散共分散行列が共通なので.
　　\(f_q(x)=\log{\dfrac{\pi_1}{\pi_2}}-\dfrac{1}{2}(x-\mu_1)'\ \Sigma_1^{-1}(x-\mu_1)+\dfrac{1}{2}(x-\mu_2)'\ \Sigma_2^{-1}(x-\mu_2)\).
即ち、
　　\(f_q(x)=\log{\dfrac{\pi_1}{\pi_2}}-\dfrac{1}{2}D_1^2+\dfrac{1}{2}D_2^2\)　　　よって、 \(f_q(x)=\log{2}+\dfrac{1}{0.75}\left[x_1+x_2\right]\). 　 これが誤判別確率最小の判別関数であり、
G1に帰属する条件は、
　　\(f_q(x)=\log{2}+\dfrac{1}{0.75}\left[x_1+x_2\right]>0\).

即ち、　\(x_1+x_2>-0.75\log{2}\).

\(x_1+x_2<-0.75\log{2}\) であれば誤判別になるので、
　　　\(c=-0.75\log2=-0.5198\)

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.81774
## [1] -0.8794717
mean(dA[,2]);mean(dB[,2])
## [1] 0.9775502
## [1] -0.9459068
#中心点
m1 <- c(mean(dA[,1]),mean(dA[,2]))
m2 <- c(mean(dB[,1]),mean(dB[,2]))
m <- (m1+m2)/2;m
## [1] -0.03086590  0.01582174
#線形判別
(Z<- lda(Z~ .,data=d))
## Call:
## lda(Z ~ ., data = d)
## 
## Prior probabilities of groups:
##   A   B 
## 0.5 0.5 
## 
## Group means:
##            X          Y
## A  0.8177400  0.9775502
## B -0.8794717 -0.9459068
## 
## Coefficients of linear discriminants:
##          LD1
## X -0.5581340
## Y -0.6061336
k <- Z $ scaling;k #法線ベクトル
##          LD1
## X -0.5581340
## Y -0.6061336
b <- m[1]*k[1]/k[2]+m[2];b #切片
## [1] -0.01259989
a <- (-1)*k[1]/k[2];a #傾き
## [1] -0.9208101
abline(b,a)

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

判別関数の値.

head(predict(Z)$x)
##          LD1
## 1 -0.4352809
## 2 -1.2781468
## 3 -1.1580889
## 4 -1.5825380
## 5 -0.7190806
## 6 -0.6596612
plot(Z,dim=2)
## Warning in rect(breaks[-n], 0, breaks[-1L], est[[grp]], col = col, ...): "dim"
## is not a graphical parameter

## Warning in rect(breaks[-n], 0, breaks[-1L], est[[grp]], col = col, ...): "dim"
## is not a graphical parameter

var(predict(Z)$x)
##          LD1
## LD1 2.116929
han<- (d[,1]+d[,2])/0.75
head(han)
## [1] 0.9693547 2.8960202 2.6745837 3.6333936 1.6631878 1.4991817
v <- (var(han[1:100])*99+var(han[101:200])*99)/198
head(han/sqrt(v))
## [1] 0.4241878 1.2672929 1.1703927 1.5899661 0.7278078 0.6560390
mean(han/sqrt(v))
## [1] -0.008777726
var(han/sqrt(v))
## [1] 2.116277
(v <- (var(han[1:100]/sqrt(v))*99+var(han[101:200]/sqrt(v))*99)/198)
## [1] 1

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

LDAの判別値は事前確率を考慮していない。

サンプル数は、200,100 とする。
hanbetu<- (d[,1]+d[,2])+0.52 par(mfrow=c(3,1)) hist(hanbetu[1:100],xlim=c(-8,8)) hist(hanbetu[101:200],xlim=c(-8,8)) hist(hanbetu,xlim=c(-8,8)) mean(hanbetu)

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(200, mA, Sigma)
dataB <- mvrnorm(100, mB, Sigma)
dA <- data.frame(X=dataA[,1],Y=dataA[,2],Z=rep("A",200))
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.04745166 0.07642428
#線形判別
(Z<- lda(Z~ .,data=d))
## Call:
## lda(Z ~ ., data = d)
## 
## Prior probabilities of groups:
##         A         B 
## 0.6666667 0.3333333 
## 
## Group means:
##            X          Y
## A  0.9712157  1.0124803
## B -0.8763124 -0.8596317
## 
## Coefficients of linear discriminants:
##          LD1
## X -0.6136808
## Y -0.4974974
k <- Z $ scaling;k
##          LD1
## X -0.6136808
## Y -0.4974974
b <- m[1]*k[1]/k[2]+m[2]+log(200/100)/k[2];b
## [1] -1.25831
b1 <- m[1]*k[1]/k[2]+m[2];b1
## [1] 0.1349576
a <- (-1)*k[1]/k[2];a
## [1] -1.233536
abline(b,a)
abline(b1,a,col="lightblue")

事前分布の差のため、判別直線は移動している。
（傾きも結構変わる印象ではある）.

判別関数の値.

head(predict(Z)$x)
##           LD1
## 1 -0.66079467
## 2 -0.08285781
## 3  0.56051918
## 4 -1.32821881
## 5 -0.06878805
## 6  0.43688315
plot(Z,dim=2)
## Warning in rect(breaks[-n], 0, breaks[-1L], est[[grp]], col = col, ...): "dim"
## is not a graphical parameter

## Warning in rect(breaks[-n], 0, breaks[-1L], est[[grp]], col = col, ...): "dim"
## is not a graphical parameter

var(predict(Z)$x)
##          LD1
## LD1 1.947581
mean(predict(Z)$x)
## [1] -1.006602e-16
han<- (d[,1]+d[,2])/0.75
head(han)
## [1]  2.39048771  1.48480793 -0.53189650  4.11985057  1.05215649 -0.01536868
var(han)
## [1] 11.25008
v <- (var(han[1:100])*99+var(han[101:200])*99)/198
mean(han/sqrt(v))
## [1] 0.4186222
var(han/sqrt(v))
## [1] 2.004434
(vp <- (var(han[1:100]/sqrt(v))*99+var(han[101:200]/sqrt(v))*99)/198)
## [1] 1
sqrt(v)
## [1] 2.369092
par(mfrow=c(2,1))
hist(-han/sqrt(v))
hist(predict(Z)$x)

var(-han/sqrt(v))
## [1] 2.004434
var(predict(Z)$x)
##          LD1
## LD1 1.947581

事後確率.

head(predict(Z)$posterior)#判別の事後確率
##           A          B
## 1 0.9409559 0.05904405
## 2 0.8285064 0.17149360
## 3 0.5612902 0.43870975
## 4 0.9844336 0.01556639
## 5 0.8243385 0.17566152
## 6 0.6228655 0.37713451
Y <- han/sqrt(v)-mean(han/sqrt(v))
head(-Y)
## [1] -0.5904089 -0.2081191  0.6431371 -1.3203775 -0.0254958  0.4251094