　「複数群の混合分布を直線で区別するのが線形判別である」
と思っていたが、そうではないようである。

データの線形一次式、例えば. $z_{1}=a_{1}x_{1}+a_{2}x_{2}$ の分散を最大とするように $\begin{pmatrix}a_1&a_2\end{pmatrix}$ を定めたとしよう。具体的には
　　　$L = w^T S_B w - \lambda (w^T S_W w - 1)$　とおいてラグランジュの未定乗数法を用いると、
　　　$S_Bw = \lambda S_W w$.
　　　即ち、$(S_W^{-1}S_B) w = \lambda w$から $w=\begin{pmatrix}a_1&a_2\end{pmatrix}:w^T S_W w =1$ が求まる。.

$z_1$ は $\begin{pmatrix}a_1&a_2\end{pmatrix}$ を方向ベクトルとする直線への正射影の長さであり.
グループが最も分離して見える直線への射影を考えたことになる.

これは主成分分析の考え方であり、重心の座標を用いて.
　　　 $z_{1}=a_{1}(x_{1}-m_1)+a_{2}(x_{2}-m_2)$ 　を判別関数とすると、
これは第１主成分に相当する。

２群の判別の場合、$z_1=0$ とおいた式が直線を表しており境界は直線になる。
この場合、 $z_1$ の正負によって所属が決定する。詳しくは、

$z_1$ の値により事後所属確率が定まり、確率の大きい群に所属が推定される。

群間分散最大化の考えは、

混合分布の分散共分散行列を群内分散と群間分散に分解する.
線形判別関数の分散は行列の２次形式なので、
分散共分散行列の分解に従い、群内分散と群間分散に分解することができる.
このとき、群間分散の割合が大きくなる射影ベクトルを定める.
射影ベクトル方向で重心を通る有向直線を判別成分とする.
これは主成分分析と同じ考えであり、第２成分まであれば２次元に投影できる.
分離をよく表す投影図において、各群の重心からの距離に従って所属を判別する.

注意.

２群の場合、射影方向を法線ベクトルとし重心を通る直線が判別直線になる.
これは判別関数を0とおいて得られる直線で判別直線と呼ばれる.
３群以上の正準判別では、各群の重心からのマハラノビス距離によって所属を判別する.
各群への事後所属確率が定まり、所属が確率的に推定される。
その際、各群の分散共分散行列が同じなら、境界は直線になる.
共通の分散共分散行列を用いるか（線形判別）、異なる分散共分散行列を用いるか（２次判別）という選択が残る.

アイスタット様「正準判別分析」の説明が分かりやすい.

Rには線形判別関数 lda (Linear Discriminant Analysis) がある。

固有ベクトルから求めた判別直線.
lda で求めた判別直線.

これらが整合していることを確かめよう.

1 2群の判別.

library(MASS) 
mA <- c(2,2)
mB <- c(-2,-2)
Sigma <- matrix(c(1,0,0,1), 2, 2)
#Sigma
dataA <- mvrnorm(100, mA, Sigma)
dataB <- mvrnorm(100, mB, matrix(c(1,0,0,1), 2, 2))
d <- rbind(dataA,dataB)
#クラス重心の中点
m1 <- c(mean(dataA[,1]),mean(dataA[,2]))
m2 <- c(mean(dataB[,1]),mean(dataB[,2]))
m <- c(mean(d[,1]),mean(d[,2]))

#群間変動
SB <- (100*(m1-m) %*% t(m1-m)+100*(m2-m) %*% t(m2-m))/200
SW <- (var(dataA)*99+var(dataB)*99)/200
SB+SW

##          [,1]     [,2]
## [1,] 5.470633 3.993211
## [2,] 3.993211 4.648252

var(d)*199/200

##          [,1]     [,2]
## [1,] 5.470633 3.993211
## [2,] 3.993211 4.648252

$S_T=S_B+S_W$ が成立していることがわかる。

dA <- data.frame(X=dataA[,1],Y=dataA[,2],Z=rep(-1,100))
dB <- data.frame(X=dataB[,1],Y=dataB[,2],Z=rep(1,100))
d <- rbind(dA,dB)
d[1,]

X	Y	Z
2.315916	2.696739	-1

2 固有値.

固有ベクトルは分散を最大とする射影方向.
射影方向は判別直線の法線ベクトル.

K <- solve(SW)  %*%  SB
eigen(K)

## eigen() decomposition
## $values
## [1] 6.765291e+00 4.089295e-16
## 
## $vectors
##           [,1]       [,2]
## [1,] 0.7084846 -0.6741270
## [2,] 0.7057263  0.7386155

k <- eigen(K)$vectors[1:2];k

## [1] 0.7084846 0.7057263

b <- m[1]*k[1]/k[2]+m[2]
a <- (-1)*k[1]/k[2]
c(a,b)#傾き、切片

## [1] -1.00390838 -0.06169791

plot(dA[,1:2],xlim=c(-5,5),ylim=c(-5,5),asp=1)
par(new=T)
plot(dB[,1:2],xlim=c(-5,5),ylim=c(-5,5),col="red",asp=1)
abline(b,a,col="blue")
abline(v = mean(d[,1]),col="lightblue")
abline(h = mean(d[,2]),col="lightblue")

2.1 判別関数の出力.

判別関数：z = a1 x1 + a2 x2 + b.

#判別関数
a <- eigen(K)$vectors[,1];a#判別係数

## [1] 0.7084846 0.7057263

#判別関数の定数項
b <- -c(mean(d[,1]),mean(d[,2]))%*%eigen(K)$vectors[,1];b

##            [,1]
## [1,] 0.04354184

a <- round(eigen(K)$vectors[,1],3);
print(判別関数 <- paste0("z = ", a[1], " * x1   ", a[2], " * x2   "," * +   ",round(b,3)))

## [1] "z = 0.708 * x1   0.706 * x2    * +   0.044"

これは、lda パッケージで求めた判別関数と若干異なっている。

3 ldaパッケージ.

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

##          LD1
## X -0.6528316
## Y -0.6502900

plot(dA[,1:2],xlim=c(-5,5),ylim=c(-5,5),asp=1)
par(new=T)
plot(dB[,1:2],xlim=c(-5,5),ylim=c(-5,5),col="red",asp=1)

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

## [1] -0.06169791

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

## [1] -1.003908

abline(b1,a1,col="red")
abline(v = mean(d[,1]),col="lightblue")
abline(h = mean(d[,2]),col="lightblue")

3.1 判別関数の出力.

判別関数：z = a1 x1 + a2 x2 + b.

#判別関数
a <- round(Z$scaling,3);a#判別係数

##      LD1
## X -0.653
## Y -0.650

Z$means

##            X         Y
## -1  2.039045  1.871019
## 1  -2.113866 -1.919301

Z$means%*%Z$scaling

##          LD1
## -1 -2.547858
## 1   2.628101

#判別関数の定数項
b <- -apply(Z$means%*%Z$scaling,2,mean);b#2列平均、判別関数の定数項

##         LD1 
## -0.04012153

print(判別関数 <- paste0("z = ", a[1], " * x1   ", a[2], " * x2   "," * +   ",round(b,3)))

## [1] "z = -0.653 * x1   -0.65 * x2    * +   -0.04"

c(mean(d[,1]),mean(d[,2]))%*%Z$scaling#ldaによる判別境界

##             LD1
## [1,] 0.04012153

c(mean(d[,1]),mean(d[,2]))%*%eigen(K)$vectors[,1]#固有値による判別境界

##             [,1]
## [1,] -0.04354184

hanbetu <- cbind(d,h =sign(a[1]*d[,1] + a[2]*d[,2] - b))
table(hanbetu$Z,hanbetu$h)

##     
##       -1   1
##   -1 100   0
##   1    0 100

3.2 結果の一覧.

(Z2<- lda(Z~ .,data=d) )

## Call:
## lda(Z ~ ., data = d)
## 
## Prior probabilities of groups:
##  -1   1 
## 0.5 0.5 
## 
## Group means:
##            X         Y
## -1  2.039045  1.871019
## 1  -2.113866 -1.919301
## 
## Coefficients of linear discriminants:
##          LD1
## X -0.6528316
## Y -0.6502900

(Z2<- lda(Z~ .,data=d,CV=T) )

## $class
##   [1] -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1
##  [26] -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1
##  [51] -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1
##  [76] -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1 -1
## [101] 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1 
## [126] 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1 
## [151] 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1 
## [176] 1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1  1 
## Levels: -1 1
## 
## $posterior
##               -1            1
## 1   1.000000e+00 3.879989e-08
## 2   9.997719e-01 2.281245e-04
## 3   9.999998e-01 1.974596e-07
## 4   1.000000e+00 3.776342e-11
## 5   1.000000e+00 2.120828e-08
## 6   1.000000e+00 2.653266e-08
## 7   9.999989e-01 1.075966e-06
## 8   9.999976e-01 2.356314e-06
## 9   9.999999e-01 9.451243e-08
## 10  9.999030e-01 9.698321e-05
## 11  9.980188e-01 1.981245e-03
## 12  9.998872e-01 1.127598e-04
## 13  9.827780e-01 1.722199e-02
## 14  9.999998e-01 1.787770e-07
## 15  9.993177e-01 6.822832e-04
## 16  9.999994e-01 5.759980e-07
## 17  9.372140e-01 6.278600e-02
## 18  1.000000e+00 8.069863e-09
## 19  9.999999e-01 7.211295e-08
## 20  9.999960e-01 3.991178e-06
## 21  9.997950e-01 2.049970e-04
## 22  9.931791e-01 6.820859e-03
## 23  1.000000e+00 9.393354e-10
## 24  9.999694e-01 3.060160e-05
## 25  9.999999e-01 1.196735e-07
## 26  1.000000e+00 8.652994e-11
## 27  9.999989e-01 1.050782e-06
## 28  9.946708e-01 5.329202e-03
## 29  1.000000e+00 1.497747e-08
## 30  1.000000e+00 5.608828e-10
## 31  9.999964e-01 3.635956e-06
## 32  9.997146e-01 2.853854e-04
## 33  9.999810e-01 1.897617e-05
## 34  9.999937e-01 6.307114e-06
## 35  1.000000e+00 2.113964e-10
## 36  9.999534e-01 4.660280e-05
## 37  1.000000e+00 2.973110e-08
## 38  9.991353e-01 8.646723e-04
## 39  9.990753e-01 9.247422e-04
## 40  9.999999e-01 7.019529e-08
## 41  9.999988e-01 1.249787e-06
## 42  9.999988e-01 1.156778e-06
## 43  9.999530e-01 4.697901e-05
## 44  9.999735e-01 2.653972e-05
## 45  9.999759e-01 2.411679e-05
## 46  9.999974e-01 2.639124e-06
## 47  9.999964e-01 3.612916e-06
## 48  9.999399e-01 6.008585e-05
## 49  9.999825e-01 1.748763e-05
## 50  9.999978e-01 2.160500e-06
## 51  9.999995e-01 4.947868e-07
## 52  9.999999e-01 1.484013e-07
## 53  9.999995e-01 4.863385e-07
## 54  1.000000e+00 2.330688e-12
## 55  9.957060e-01 4.293968e-03
## 56  9.999931e-01 6.889787e-06
## 57  9.999031e-01 9.687782e-05
## 58  9.997167e-01 2.833014e-04
## 59  1.000000e+00 6.730200e-09
## 60  9.999887e-01 1.130582e-05
## 61  9.999911e-01 8.852755e-06
## 62  9.999928e-01 7.165941e-06
## 63  9.999938e-01 6.213584e-06
## 64  9.999999e-01 9.864420e-08
## 65  1.000000e+00 1.161311e-09
## 66  9.999172e-01 8.280637e-05
## 67  1.000000e+00 8.013367e-09
## 68  9.999992e-01 8.479259e-07
## 69  1.000000e+00 4.223377e-08
## 70  9.999567e-01 4.333192e-05
## 71  9.999992e-01 7.848887e-07
## 72  9.982981e-01 1.701887e-03
## 73  9.999998e-01 2.333800e-07
## 74  1.000000e+00 4.152036e-10
## 75  1.000000e+00 1.630195e-11
## 76  9.786011e-01 2.139890e-02
## 77  9.999279e-01 7.210292e-05
## 78  9.999513e-01 4.868956e-05
## 79  9.999999e-01 9.743812e-08
## 80  9.999680e-01 3.195230e-05
## 81  9.997983e-01 2.016693e-04
## 82  9.997997e-01 2.002709e-04
## 83  9.999977e-01 2.337426e-06
## 84  9.999989e-01 1.065182e-06
## 85  1.000000e+00 5.660170e-09
## 86  9.999986e-01 1.411635e-06
## 87  9.981057e-01 1.894261e-03
## 88  1.000000e+00 1.268488e-08
## 89  9.999865e-01 1.349214e-05
## 90  9.999850e-01 1.497658e-05
## 91  1.000000e+00 1.845414e-10
## 92  1.000000e+00 2.934118e-10
## 93  9.999858e-01 1.415069e-05
## 94  9.999986e-01 1.401849e-06
## 95  1.000000e+00 2.339150e-09
## 96  9.999796e-01 2.037425e-05
## 97  1.000000e+00 2.664222e-09
## 98  9.812882e-01 1.871181e-02
## 99  9.999997e-01 2.640614e-07
## 100 1.000000e+00 1.289056e-11
## 101 1.274952e-03 9.987250e-01
## 102 5.802552e-08 9.999999e-01
## 103 7.737773e-06 9.999923e-01
## 104 3.965515e-08 1.000000e+00
## 105 7.384466e-06 9.999926e-01
## 106 5.886928e-06 9.999941e-01
## 107 4.292999e-06 9.999957e-01
## 108 2.128264e-07 9.999998e-01
## 109 3.334388e-06 9.999967e-01
## 110 2.949814e-08 1.000000e+00
## 111 1.353452e-08 1.000000e+00
## 112 1.721813e-03 9.982782e-01
## 113 3.714253e-07 9.999996e-01
## 114 2.952939e-07 9.999997e-01
## 115 1.558233e-05 9.999844e-01
## 116 2.093214e-12 1.000000e+00
## 117 5.641511e-08 9.999999e-01
## 118 2.025858e-09 1.000000e+00
## 119 7.200938e-09 1.000000e+00
## 120 5.248601e-03 9.947514e-01
## 121 6.927461e-08 9.999999e-01
## 122 4.991685e-06 9.999950e-01
## 123 7.152522e-06 9.999928e-01
## 124 1.267895e-05 9.999873e-01
## 125 1.167747e-07 9.999999e-01
## 126 1.847508e-07 9.999998e-01
## 127 6.921311e-04 9.993079e-01
## 128 2.551286e-02 9.744871e-01
## 129 1.366860e-06 9.999986e-01
## 130 4.186198e-08 1.000000e+00
## 131 1.824295e-05 9.999818e-01
## 132 3.683172e-06 9.999963e-01
## 133 2.710319e-06 9.999973e-01
## 134 9.749963e-09 1.000000e+00
## 135 9.519731e-06 9.999905e-01
## 136 3.968711e-07 9.999996e-01
## 137 2.385130e-01 7.614870e-01
## 138 4.301493e-07 9.999996e-01
## 139 5.563126e-05 9.999444e-01
## 140 1.180460e-03 9.988195e-01
## 141 2.088407e-06 9.999979e-01
## 142 2.846805e-09 1.000000e+00
## 143 5.115445e-04 9.994885e-01
## 144 3.632888e-05 9.999637e-01
## 145 2.124153e-06 9.999979e-01
## 146 2.060231e-06 9.999979e-01
## 147 9.193325e-06 9.999908e-01
## 148 5.927436e-05 9.999407e-01
## 149 1.255956e-10 1.000000e+00
## 150 7.788567e-07 9.999992e-01
## 151 1.252929e-04 9.998747e-01
## 152 1.003837e-04 9.998996e-01
## 153 1.253482e-06 9.999987e-01
## 154 3.473820e-04 9.996526e-01
## 155 8.134541e-09 1.000000e+00
## 156 3.986992e-04 9.996013e-01
## 157 3.496601e-08 1.000000e+00
## 158 1.235656e-06 9.999988e-01
## 159 1.357482e-08 1.000000e+00
## 160 4.373278e-05 9.999563e-01
## 161 1.485635e-09 1.000000e+00
## 162 5.898504e-07 9.999994e-01
## 163 1.964166e-07 9.999998e-01
## 164 1.296503e-04 9.998703e-01
## 165 2.297669e-04 9.997702e-01
## 166 9.068682e-08 9.999999e-01
## 167 2.292193e-05 9.999771e-01
## 168 3.166783e-07 9.999997e-01
## 169 2.443203e-07 9.999998e-01
## 170 7.002415e-09 1.000000e+00
## 171 7.132763e-04 9.992867e-01
## 172 2.761928e-07 9.999997e-01
## 173 6.548296e-05 9.999345e-01
## 174 1.981953e-06 9.999980e-01
## 175 1.334573e-03 9.986654e-01
## 176 1.809414e-03 9.981906e-01
## 177 6.844892e-11 1.000000e+00
## 178 4.457976e-11 1.000000e+00
## 179 5.118228e-12 1.000000e+00
## 180 2.619995e-06 9.999974e-01
## 181 1.827740e-04 9.998172e-01
## 182 4.743575e-09 1.000000e+00
## 183 6.586822e-07 9.999993e-01
## 184 6.095779e-04 9.993904e-01
## 185 5.245548e-07 9.999995e-01
## 186 5.608125e-09 1.000000e+00
## 187 1.362509e-08 1.000000e+00
## 188 6.591882e-06 9.999934e-01
## 189 1.338416e-08 1.000000e+00
## 190 6.491232e-08 9.999999e-01
## 191 8.273443e-04 9.991727e-01
## 192 8.441910e-04 9.991558e-01
## 193 9.276367e-08 9.999999e-01
## 194 6.012388e-05 9.999399e-01
## 195 3.565619e-04 9.996434e-01
## 196 4.625005e-02 9.537500e-01
## 197 5.072336e-11 1.000000e+00
## 198 3.917158e-02 9.608284e-01
## 199 2.229799e-11 1.000000e+00
## 200 1.769011e-07 9.999998e-01
## 
## $terms
## Z ~ X + Y
## attr(,"variables")
## list(Z, X, Y)
## attr(,"factors")
##   X Y
## Z 0 0
## X 1 0
## Y 0 1
## attr(,"term.labels")
## [1] "X" "Y"
## attr(,"order")
## [1] 1 1
## attr(,"intercept")
## [1] 1
## attr(,"response")
## [1] 1
## attr(,".Environment")
## <environment: R_GlobalEnv>
## attr(,"predvars")
## list(Z, X, Y)
## attr(,"dataClasses")
##         Z         X         Y 
## "numeric" "numeric" "numeric" 
## 
## $call
## lda(formula = Z ~ ., data = d, CV = T)
## 
## $xlevels
## named list()

d$判定 <- a[1]*d$X + a[2]*d$Y - b 

d$エラー <- ifelse(d$Z * d$判定 >= 0, "　", "×")
d["正群確率"] <- Z2$posterior[,2]
d["負群確率"] <- Z2$posterior[,1]
head(d)

X	Y	Z	判定	正群確率	負群確率
2.3159159	2.696739	-1	-3.225052	0.0000000	1.0000000
0.6125035	1.826977	-1	-1.547378	0.0002281	0.9997719
3.5052257	1.025314	-1	-2.915245	0.0000002	0.9999998
4.2304230	2.717113	-1	-4.488468	0.0000000	1.0000000
2.3869783	2.800825	-1	-3.339111	0.0000000	1.0000000
3.3693510	1.747321	-1	-3.295823	0.0000000	1.0000000

#Z2

hist(d$判定)

head(d$判定)

## [1] -3.225052 -1.547378 -2.915245 -4.488468 -3.339111 -3.295823

mean(d$判定)

## [1] 0.08024237

sd(d$判定[1:100])

## [1] 1.003139

sd(d$判定[101:200])#群内分散１

## [1] 0.9966958

p1 <- dnorm(predict(Z)$x[1:10], mean=mean(d$判定[1:100]), sd=1);p1

##  [1] 0.29012013 0.27098570 0.35421764 0.04767667 0.26315343 0.27360848
##  [7] 0.39355977 0.39885157 0.32779859 0.31253357

p2 <- dnorm(predict(Z)$x[1:10], mean=mean(d$判定[101:200]), sd=1);p2

##  [1] 7.111595e-09 3.919163e-05 4.330689e-08 1.692798e-12 3.574243e-09
##  [6] 4.660827e-09 2.558191e-07 5.448557e-07 1.950834e-08 1.773478e-05

p1/(p1+p2)

##  [1] 1.0000000 0.9998554 0.9999999 1.0000000 1.0000000 1.0000000 0.9999993
##  [8] 0.9999986 0.9999999 0.9999433

p2/(p1+p2)

##  [1] 2.451259e-08 1.446053e-04 1.222607e-07 3.550580e-11 1.358235e-08
##  [6] 1.703466e-08 6.500129e-07 1.366059e-06 5.951318e-08 5.674197e-05

#str(Z2)

固有ベクトルで求めた境界と lda で求めた境界は、ほぼ一致している。
少し異なっているが、理由は分からない。

個数を変える.
分散を変える.

データの判別

4 ３群の正準判別.

３群の線形判別は、正準判別と呼ばれる。
大森宏先生の説明によると、

u1 が最もよく群間の変動を説明し，
u2 がその次に群間の変動を説明する，

というようになる．この ui を正準変量と呼ぶこともある．
　通常は，u1，u2 で群間の変動のかなりの部分が説明されるので，データをこの２つの正順得点で (u1，u2 ) 平面上にプロットし，判別関数で平面を分割すれば，判別の様子が視覚化される．
　また，u1 と u2 は群内で無相関となる．しかし， ai’aj ≠ 0 なので，ui 軸と uj軸は直交しない．この点が主成分分析による次元の減少とは異なっている.

正準得点による群内平方和積和行列は単位行列になるので，未分類の対象データの正準得点 u ＝ (u1，u2)’ が与えられたとすると，グループ j の平均とのマハラノビス距離は，

　　　$D^2=(u-\bar{u})'(u-\bar{u})$

のようになり，(U1，U2) 平面上では Euclid 距離となる．

大森先生の説明によると、３群の正準判別では.
(u1，u2 )のマハラノビス距離（この場合はユークリッド距離になっている）に従って、３本の線形判別関数が求まる。

ここでは主成分方向を図示するにとどめる。

library(MASS) 
mC <- c(-2,2)
dataC <- mvrnorm(100, mC, Sigma)
d <- rbind(dataA,dataB,dataC)
#var(dataA)*99/100
#var(dataB)*99/100
#var(d)*299/300
#クラス重心の中点
m3 <- c(mean(dataC[,1]),mean(dataC[,2]))
m <- c(mean(d[,1]),mean(d[,2]))
#m
#群間変動
S1 <- (m1-m) %*% t(m1-m)
S2 <- (m2-m) %*% t(m2-m)
S3 <- (m3-m) %*% t(m3-m)
SB <- (100*S1+100*S2+100*S3)/300
SW <- (var(dataA)*99+var(dataB)*99+var(dataC)*99)/300
SB+SW

##          [,1]     [,2]
## [1,] 4.762756 1.774136
## [2,] 1.774136 4.482228

var(d)*299/300

##          [,1]     [,2]
## [1,] 4.762756 1.774136
## [2,] 1.774136 4.482228

$S_T=S_B+S_W$ が成立していることがわかる。

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))
dC <- data.frame(X=dataC[,1],Y=dataC[,2],Z=rep("C",100))
d <- rbind(dA,dB,dC)

4.1 固有値.

K <- solve(SW)  %*%  SB
eigen(K)

## eigen() decomposition
## $values
## [1] 4.652881 1.848766
## 
## $vectors
##           [,1]       [,2]
## [1,] 0.7623175 -0.6648356
## [2,] 0.6472032  0.7469897

k1 <- eigen(K)$vectors[1:2];k1

## [1] 0.7623175 0.6472032

k2 <- eigen(K)$vectors[3:4];k2

## [1] -0.6648356  0.7469897

plot(dA[,1:2],xlim=c(-5,5),ylim=c(-5,5),asp=1)
par(new=T)
plot(dB[,1:2],xlim=c(-5,5),ylim=c(-5,5),col="red",asp=1)
par(new=T)
plot(dC[,1:2],xlim=c(-5,5),ylim=c(-5,5),col="blue",asp=1)
b1 <- m[1]*k1[1]/k1[2]+m[2];b

##         LD1 
## -0.04012153

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

##      LD1
## X -0.653
## Y -0.650

abline(b1,a1)

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

##         LD1 
## -0.04012153

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

##      LD1
## X -0.653
## Y -0.650

abline(b2,a2,col="blue")

4.2 LDA.

(Z<- lda(Z~ .,data=d))

## Call:
## lda(Z ~ ., data = d)
## 
## Prior probabilities of groups:
##         A         B         C 
## 0.3333333 0.3333333 0.3333333 
## 
## Group means:
##           X         Y
## A  2.039045  1.871019
## B -2.113866 -1.919301
## C -1.979575  2.104534
## 
## Coefficients of linear discriminants:
##          LD1        LD2
## X -0.7130802 -0.6620500
## Y -0.6054009  0.7438599
## 
## Proportion of trace:
##    LD1    LD2 
## 0.7156 0.2844

#Z
plot(Z,dimen=2)

k <- Z $ scaling;k

##          LD1        LD2
## X -0.7130802 -0.6620500
## Y -0.6054009  0.7438599

apply(Z$means %*% Z$scaling, 2, mean)

##        LD1        LD2 
## 0.07336409 0.96322532

k1 <- eigen(K)$vectors[1:2];k1

## [1] 0.7623175 0.6472032

k2 <- eigen(K)$vectors[3:4];k2

## [1] -0.6648356  0.7469897

plot(dA[,1:2],xlim=c(-5,5),ylim=c(-5,5),asp=1)
par(new=T)
plot(dB[,1:2],xlim=c(-5,5),ylim=c(-5,5),col="red",asp=1)
par(new=T)
plot(dC[,1:2],xlim=c(-5,5),ylim=c(-5,5),col="blue",asp=1)
b1 <- m[1]*k1[1]/k1[2]+m[2];b

##         LD1 
## -0.04012153

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

##      LD1
## X -0.653
## Y -0.650

abline(b1,a1,col="red")

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

##         LD1 
## -0.04012153

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

##      LD1
## X -0.653
## Y -0.650

abline(b2,a2,col="red")

判別結果.

①　 predict(Z)$class.
②　 predict(Z)$posterior.
③　 predict(Z)$x.

predict(Z)$class

##   [1] A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A A
##  [38] A A A A A A A A A A A A A A A A A A A A A A A A A C A A A A A A A A A A A
##  [75] A A A A A A A A A A A A A A A A A A A A A A A A A A B B B B B B B B B B B
## [112] B B B B B B B B C B B B B B B B B B B B B B B B B B B B B B B B B B B B B
## [149] B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B B C B
## [186] B B B B B B C B B B B B B B B C B C C C C C C C C C C C C C C C C C C C C
## [223] C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C C
## [260] C C C C C C C C C C C C C C C C C C C C C C C B C C C C C C C C C C C C C
## [297] C C C C
## Levels: A B C

head(predict(Z)$posterior)

##           A            B            C
## 1 0.9997764 2.483820e-08 2.235873e-04
## 2 0.9086402 2.058794e-04 9.115388e-02
## 3 0.9999988 6.774854e-08 1.106098e-06
## 4 0.9999998 1.944776e-11 1.507525e-07
## 5 0.9998212 1.363664e-08 1.787945e-04
## 6 0.9999974 1.084515e-08 2.588704e-06

head(predict(Z)$x)

##         LD1          LD2
## 1 -3.357406 -0.490481161
## 2 -1.616182 -0.009718415
## 3 -3.193597 -2.521169532
## 4 -4.734937 -1.742825558
## 5 -3.471093 -0.460103029
## 6 -3.533811 -1.894142309

4.3 群内分散.

library("dplyr")

## 
## Attaching package: 'dplyr'

## 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	LD2
-3.357406	-0.4904812
-1.616182	-0.0097184
-3.193597	-2.5211695
-4.734937	-1.7428256
-3.471093	-0.4601030
-3.533811	-1.8941423

dat <- mutate(predict,Z=c(rep(-1,100),rep(0,100),rep(1,100)))
head(dat)

LD1	LD2	Z
-3.357406	-0.4904812	-1
-1.616182	-0.0097184	-1
-3.193597	-2.5211695	-1
-4.734937	-1.7428256	-1
-3.471093	-0.4601030	-1
-3.533811	-1.8941423	-1

Z<- lda(Z~ .,data=dat)
plot(Z,dimen=2)
abline(0,1)

var(predict[1:100,])

##            LD1        LD2
## LD1 1.03896325 0.03696072
## LD2 0.03696072 1.15883217

var(predict[101:200,])

##            LD1        LD2
## LD1 1.03507607 0.09191274
## LD2 0.09191274 0.93313096

var(predict[201:300,])

##            LD1        LD2
## LD1  0.9259607 -0.1288735
## LD2 -0.1288735  0.9080369

#プール分散
(var(predict[1:100,])*99+var(predict[101:200,])*99+var(predict[201:300,])*99)/300

##               LD1           LD2
## LD1  9.900000e-01 -2.782959e-16
## LD2 -2.782959e-16  9.900000e-01

プールされた群内分散は単位行列になっている。
mahalanobis(dat,m.dat,v.dat).

mA <- apply(predict[1:100,],2,mean);mA

##        LD1        LD2 
## -2.6600834 -0.9213991

mB <- apply(predict[101:200,],2,mean);mB

##        LD1        LD2 
##  2.5959388 -0.9914318

mC <- apply(predict[201:300,],2,mean);mC

##        LD1        LD2 
## 0.06414462 1.91283099

p1 <- mahalanobis(predict[1:6,],mA,diag(c(1,1)))
P1 <- exp(-1/2*p1);P1

##          1          2          3          4          5          6 
## 0.71464116 0.38272207 0.24124284 0.08291989 0.64709298 0.42536202

p2 <- mahalanobis(predict[1:6,],mB,diag(c(1,1)))
P2 <- exp(-1/2*p2);P2

##            1            2            3            4            5            6 
## 1.775437e-08 8.671704e-05 1.634387e-08 1.612607e-12 8.825753e-09 4.613128e-09

p3 <- mahalanobis(predict[1:6,],mC,diag(c(1,1)))
P3 <- exp(-1/2*p3);P3

##            1            2            3            4            5            6 
## 1.598204e-04 3.839429e-02 2.668386e-07 1.250038e-08 1.157173e-04 1.101139e-06

P1/(P1+P2+P3)

##         1         2         3         4         5         6 
## 0.9997764 0.9086402 0.9999988 0.9999998 0.9998212 0.9999974

4.4 まとめ.

線形判別分析は、
　各群の正規性と等分散性を仮定している.

２群の場合.
線形判別関数によって１次元に投影して判別する。
投影方向は任意であるが、分散が最大の方向が識別が良いと考えられる.
それは主成分の方向ということになる.
等分散を仮定しているため、
　　　群内分散はプール分散を採用する.
判別関数LD1の値は混合分布を形成する.
　　　各群でのLD1の群内分散は１に近く.
　　　プールした分散は１になる.
　　　プール分散によるマハラノビスの距離から帰属確率が定まり、
事後確率は、

$\begin{split} P(A|x)&=\dfrac{\pi_Af_A(x)}{\pi_Af_A(x)+\pi_Bf_B(x)} \end{split}$.

などによってが決まる.

３群の場合.
線形判別関数によって１次元に投影して判別する。
分散の大きさに応じて、2つの投影方向が定まる.
それは主成分の方向である.
(LD1,LD2)の２次元混合分布について、
　　　プール分散は単位行列である.
　　　プール分散による各群の重心とのマハラノビスの距離から帰属確率が定まる.
事後確率は、
　　　$\begin{split} P(A|x)&=\dfrac{\pi_Af_A(x)}{\pi_Af_A(x)+\pi_Bf_B(x)+\pi_Cf_C(x)} \end{split}$.

などによって決まる（上記の例を参照）.
これは、RのLDRで求めた事後確率と一致する。

5 人文2018③.

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)
d <- rbind(dataA,dataB)
var(dataA)*99

##           [,1]      [,2]
## [1,] 1638.6682  910.1252
## [2,]  910.1252 1327.9505

var(dataB)*99

##          [,1]     [,2]
## [1,] 1905.822 1207.936
## [2,] 1207.936 1894.243

var(d)*199

##          [,1]     [,2]
## [1,] 7879.736 1670.817
## [2,] 1670.817 3268.333

#クラス重心の中点
m1 <- c(mean(dataA[,1]),mean(dataA[,2]))
m2 <- c(mean(dataB[,1]),mean(dataB[,2]))
m <- c(mean(d[,1]),mean(d[,2]))
m

## [1] 34.68246 49.99685

#群間変動
print((m1-m) %*% t(m1-m))

##           [,1]       [,2]
## [1,] 21.676229 -2.2362242
## [2,] -2.236224  0.2306997

print((m2-m) %*% t(m2-m))

##           [,1]       [,2]
## [1,] 21.676229 -2.2362242
## [2,] -2.236224  0.2306997

SB <- 99*(m1-m) %*% t(m1-m)+99*(m2-m) %*% t(m2-m)
SW <- var(dataA)*99+var(dataB)*99
SB+SW

##          [,1]     [,2]
## [1,] 7836.383 1675.289
## [2,] 1675.289 3267.872

var(d)*199

##          [,1]     [,2]
## [1,] 7879.736 1670.817
## [2,] 1670.817 3268.333

(SB/100+SW/100)*0.5

##           [,1]      [,2]
## [1,] 39.181916  8.376445
## [2,]  8.376445 16.339358

var(d)*199/200

##           [,1]      [,2]
## [1,] 39.398679  8.354083
## [2,]  8.354083 16.341665

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)

head(dA)

X	Y	Z
30.44660	48.43592	A
35.55513	56.73468	A
33.20638	56.83695	A
34.12411	56.13370	A
32.75468	51.46630	A
22.12735	48.36094	A

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

## [1] 34.68246 49.99685

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

##          LD1
## X  0.3023711
## Y -0.2182716

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

## [1] 1.95133

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

## [1] 1.385298

abline(b,a)

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.21960 50.44154

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

##          LD1
## X  0.3234955
## Y -0.2177155

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

## [1] -1.889979

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

## [1] 1.485864

abline(b,a)

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

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

## [1] 1.397157

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

5.0.2 ldaの結果.

## Call:
## lda(Z ~ ., data = d)
## 
## Prior probabilities of groups:
##   A   B 
## 0.5 0.5 
## 
## Group means:
##          X        Y
## A 30.63988 50.58600
## B 39.79932 50.29709
## 
## Coefficients of linear discriminants:
##          LD1
## X  0.3234955
## Y -0.2177155

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

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

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

##      A    B
## 1 0.14 0.86
## 2 1.00 0.00
## 3 0.82 0.18
## 4 0.99 0.01
## 5 0.96 0.04
## 6 0.98 0.02

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

##          LD1
## 1  0.5871917
## 2 -2.3364358
## 3 -0.4907372
## 4 -1.3921542
## 5 -1.0140406
## 6 -1.3590450