多次元の判別分析について整理する。
判別分析は混合分布に対し、その所属を判別する。
各群の重心からの確率距離を考え、事後確率に従って判別することになる。
まとめとしては、東京大学大森宏先生のレジメが最も良い。
大森先生の資料に従ってまとめてみよう。
大森宏「判別分析」.

1 判別分析とは.

データを判別する手法である．分類のわかっているデータ（トレイニングデータ，教師データ）があり，データの属性値に基づいて，この分類を判別する手法である．機械学習の分野では，教師あり学習（supervised learning）と呼ばれている．データ分類手法として知られるクラスター分析は，分類がわかっていないので，教師無し学習（unsupervised learning）と呼ばれる．
　応用分野は多岐にわたり，病気診断，スニップデータからの発現分類，スパムメールフィルター，などがあげられる．

2 ２群の判別.

2.1 空間の分割と判別規則.

母集団が２つの部分母集団 A，B の混合からなっている場合を考える．部分母集団 A，B の確率密度をそれぞれ fA(x )，fB(x ) とし，その混合比率をπA，πB とすると，母集団全体の確率密度は，
　　　f(x ) ＝ πAfA(x ) ＋ πBfB(x ).
となる．いま，帰属未知の対象 X の属性値 x に基づき部分母集団への帰属を決めるとする．
このとき， x の取り得る値の全体（標本空間） R　は，RA と RB に分割され，
　　x ∈ RA なら X は A に帰属，　x ∈ RB なら X は B に帰属．
となる．これを判別規則という．

2.2 最適な判別規則.

　x の空間 R の分割が判別規則に対応するが，どのような空間分割を行うのがよいであろうか．まず，ベイズの定理による事後確率の最大化という考え方がある．

ベイズ（Bayes）の定理.
P(A) を事象 A が起こる確率，P(A|B) を事象 B が生起したもとでの事象 A が起こる条件付き確率， P(A∩B) を事象 A と事象 B が同時に生起する確率とすると，
　　　P(A∩B) ＝ P(A|B)P(B) ＝ P(B|A)P(A).
が成り立つので.
　　　P(A|B) ＝ P(B|A)P(A)/P(B)
である．これをベイズの定理と言う．
いま，標本空間 Ω が互いに排反な事象，
　　　A1，…，An 　（P(Ai∩Aj，i≠j， P(A1∪…∪An) ＝ P(A1)＋…＋P(An)＝1），
に分割されたとすると，Ω の任意の部分集合 B の生起する確率 P(B) は，
　　　P(B) ＝ P(B|A1)P(A1)＋…＋P(B|An)P(An) ＝ ΣiP(B|Ai)P(Ai) と表現される．
これより，事象 B が生起したときの事象 Ai が生起する事後確率は，
　　　\(\begin{split} P(A_i|B)&=\dfrac{P(B|A_i)P(A_i)}{P(B)}=\dfrac{P(B|A_i)P(A_i)}{\sum_jP(B|A_j)P(A_j)} \end{split}\).
と表せる．

ベイズ判別規則.
　対象 X が A に帰属したとき，対象 X の値が x のごく近傍を取る確率は， P(x |A) ＝ fA(x )Δx であり，対象 X が B に帰属したとき， x のごく近傍を取る確率は，P(x |B) ＝ fB(x )Δx とあらわせる．　対象 X の値が x の近傍であるとき，X が部分母集団 A 由来である事後確率は， P(A) ＝ πA に注意すると，

　　　\(\begin{split} P(A|x)&=\dfrac{P(x|A)P(A)}{P(x)}=\dfrac{f_A(x)\Delta x\pi_A}{P(x|A)P(A)+P(x|B)P(B)}\\ &=\dfrac{f_A(x)\Delta x\pi_A}{f_A(x)\Delta x\pi_A+f_B(x)\Delta x\pi_B}=\dfrac{\pi_Af_A(x)}{\pi_Af_A(x)+\pi_Bf_B(x)} \end{split}\).
であり，X が部分母集団 B 由来である事後確率は，
　　　\(\begin{split} P(B|x)&=\dfrac{\pi_Bf_B(x)}{\pi_Af_A(x)+\pi_Bf_B(x)} \end{split}\).
となる．事後確率の高い部分集団に対象 X を帰属させるとすると，

\(X\) が部分母集団 Aに由来と判定する条件は、
　　　\(P(A|x)>P(B|x) \Leftrightarrow \pi_Af_A(x)>\pi_Bf_B(x)\).

　　\(\therefore R_A=\{x|\pi_Af_A(x)>\pi_Bf_B(x)\}\).
即ち，
　　　RA ＝ { x | logπA ＋ logfA(x ) ＞ logπB ＋ logfB(x ) }

判別関数.
　　　\(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\)　と表される。.

誤判別確率最小判別.
判別規則により，対象Xを判別したとする．このときは次のような誤判別がある．

X は A から得られたのに B から得られたと判別する．この確率を Pr[B|A] とおく．
逆に，X は B から得られたのに A から得られたと判別する．この確率を Pr[A|B] とおく．

すると，誤判別の確率は，

　　　\(\begin{split} 　　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\\ &=\pi_A\left[\int_{R}f_A(x)dx-\int_{R_A}f_A(x)dx\right]+\pi_B\int_{R_A}f_B(x)dx\\ &=\pi_A+\int_{R_A}\left[\pi_Bf_B(x)-\pi_Af_A(x)\right]dx \end{split}\).
となる．これより，RA として，
　　　\(\pi_Bf_B(x)-\pi_Af_A(x)<0\).　の部分のみを含み，
　　　\(\pi_Bf_B(x)-\pi_Af_A(x)>0\).　の部分をふくまなければ.

誤判別確率は最小となる．

次のように説明してもよい。

\(T(R_A,R_B)=\dfrac{1}{2}+\dfrac{1}{2}\int_{R_A}\left[f_B(x)-f_A(x)\right]dx\) が最小となるのは、
　　　　\(R_A\)において \(f_A(x)\geq f_B(x)\) のとき.
\(T(R_A,R_B)=\dfrac{1}{2}+\dfrac{1}{2}\int_{R_B}\left[f_A(x)-f_B(x)\right]dx\) が最小となるのは、
　　　　\(R_B\)において \(f_A(x)\leq f_B(x)\) のとき.

これを満たすのは、\(R_A,R_B\)の境界が、\(f_A(x),f_B(x)\) の交点の射影する点であるときである。
（２つの群の分布円の交線によって境界がわかたれるときである？）.
人文2018②参照.

相関比 correlation ratio.
2変数のうち一方が量的変数、もう一方が名義尺度である場合に、この2つの変数間の相関の程度を表す統計量。
0から1までの値をとる.
統計学がわかった様「相関比の意味と計算方法」参照.
（質的データを扱う場合は、相関係数の計算はできません）.

名義尺度に対する変動が大きければ相関比は大きく.
名義尺度に対する変動が小さければ相関比は小さい.

従って、相関比を最大化すれば群間変動は最大化され、判別しやすくなる.

相関比.

2.3 正規母集団のもとでの判別.

２次判別.
部分母集団が正規分布で記述されているときは，簡単な判別規則をつくることができる．
平均 μ，分散共分散行列 Σ の２次元正規分布の密度関数は，
　　　\(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]\).

　　　\(\Sigma=\begin{pmatrix} \sigma_{x}^{2} & \sigma_{xy}\\ \sigma_{xy} & \sigma_{y}^{2} \end{pmatrix}\).

　　　\(D^2=\begin{pmatrix}x-\mu_{x}&y-\mu_{y}\end{pmatrix}\begin{pmatrix}\sigma_{x}^{2} &\sigma_{xy}\\ \sigma_{xy} &\sigma_{y}^{2}\end{pmatrix}^{-1} \begin{pmatrix}x-\mu_x\\y-\mu_y\end{pmatrix}\).
　

準１級2019_論述③ 参照（_直前問題_判別分析、類題　ワーク問23_3）.

２次判別における判別器.

　sign(符号関数)：実数に対しその符号に応じて1、−1、0のいずれかを返す関数。

\(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\).

　　　\(f_q(x)>0 \Leftrightarrow D_1^2 < D_2^2\).

マハラノビスの距離の近い群に所属することになる。
　

さて、大森先生のレジメに戻ろう。

A,B が正規母集団で \(\pi_A=\pi_B=\dfrac{1}{2}\)　かつ分散共分散行列が等しければ、
　　　\(P(A|x)>P(B|x) \Leftrightarrow \exp\left({-\dfrac{1}{2}}D_A^2\right)>\exp\left({-\dfrac{1}{2}}D_B^2\right) \Leftrightarrow D_A^2 < D_B^2\).
　　　\(\therefore R_A=\{x|D_A^2 < D_B^2\}\).

判別器　\(f_q(x)\)を用いれば、
　　　\(f_q(x)=-\dfrac{1}{2}D_A^2+\dfrac{1}{2}D_B^2=\dfrac{1}{2}(D_B^2-D_A^2)\).

これは大森先生のレジメの 線形判別関数（linear discriminant function） \(w\) と一致している。

　　　\(D_B^2-D_A^2=(x-\mu_B)'\Sigma^{-1}(x-\mu_B)-(x-\mu_A)'\Sigma^{-1}(x-\mu_A)\).

展開は \(\Sigma^{-1}\) を間に挟むだけで、通常の展開と同型になる。

　　　\(\begin{split} (x-\mu_B)'(x-\mu_B)-(x-\mu_A)'(x-\mu_A)&= x'x-x'\mu_B-\mu_B'x+\mu_B'\mu_B-x'x+x'\mu_A+\mu_A'x-\mu_A'\mu_A\\ &=-x'\mu_B-\mu_B'x+\mu_B'\mu_B+x'\mu_A+\mu_A'x-\mu_A'\mu_A\\ &=-x'(\mu_B-\mu_A)-(\mu_B-\mu_A)'x+\mu_B'\mu_B-\mu_A'\mu_A\\ &=-2x'(\mu_B-\mu_A)+(\mu_B'+\mu_A')(\mu_B-\mu_A) \end{split}\).

　　　\(\begin{split} D_B^2-D_A^2&=-2\cdot x'\Sigma^{-1}(\mu_B-\mu_A)+(\mu_B'+\mu_A')\Sigma^{-1}(\mu_B-\mu_A)\\ &=-2\left(x-\dfrac{\mu_B+\mu_A}{2}\right)'\Sigma^{-1}(\mu_B-\mu_A) \end{split}\).

よって、
　　　\(w=f_q(x)=\left(x-\dfrac{\mu_A+\mu_B}{2}\right)'\Sigma^{-1}(\mu_A-\mu_B)\)

\(D_B^2-D_A^2=0\) とすると、

　　　\(\left(x-\dfrac{\mu_B+\mu_A}{2}\right)'\Sigma^{-1}(\mu_B-\mu_A)=0\).

これは、グループA,Bの平均の中点を通り、\(\Sigma^{-1}(\mu_B-\mu_A)\) を法線ベクトルとする直線になる。

直線の法線ベクトルを \(\begin{pmatrix} a\\b \end{pmatrix}\) とおくと、 \(\begin{pmatrix} a\\b \end{pmatrix}=\Sigma^{-1}(\mu_B-\mu_A)\)　　\(\therefore \Sigma\begin{pmatrix} a\\b \end{pmatrix}=(\mu_B-\mu_A)\).

分散共分散行列を係数とし平均の差を定数とする連立方程式を解くと \(a,b\) が求まる。

2018人文② 参照.

3 線形判別の例.

3.1 線形判別関数の例１.

ひょっこり統計線形判別を参考.

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

大森宏「判別分析」の例による.
サンプル数を100ずつとする。

library(MASS) 
mA <- c(50,30)
mB <- c(70,24)
Sigma <- matrix(c(100,35,35,25), 2, 2)
Sigma

##      [,1] [,2]
## [1,]  100   35
## [2,]   35   25

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,100),ylim=c(0,50))
par(new=T)
plot(dB[,1:2],xlim=c(20,100),ylim=c(0,50),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] 59.65979 26.96378

par(new=T)
plot(x=mean(dA[,1]),y=mean(dA[,2]),xlim=c(20,100),ylim=c(0,50),pch=15,col="blue",xlab="",ylab="")
par(new=T)
plot(x=mean(dB[,1]),y=mean(dB[,2]),xlim=c(20,100),ylim=c(0,50),pch=15,col="blue",xlab="",ylab="")
par(new=T)
plot(m[1],m[2],xlim=c(20,100),ylim=c(0,50),pch=17,col="blue",xlab="",ylab="")

(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 49.93114 29.82596
## B 69.38844 24.10160
## 
## Coefficients of linear discriminants:
##          LD1
## X  0.1364599
## Y -0.2488984

# ldaから求めた法線ベクトル
k <- Z $ scaling;norm(k)

## [1] 0.3853584

k/norm(k)

##          LD1
## X  0.3541118
## Y -0.6458882

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

## [1] -5.745025

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

## [1] 0.5482555

abline(b,a)

3.1.1.1 参考.

ax+by は.

w=(a,b)のノルムが１のとき、(a,b) 方向への (x,y) の正射影の長さを意味する.
w=(a,b)のノルムが一般のとき、(a,b) 方向への (x,y) の正射影の長さに比例する.

線形判別では、
　\(w^T S_W w=1\)という条件のもとで.
　　　\({J(w) = \dfrac{w^T S_B w}{w^T S_W w} }=w^T S_B w\) の最大を考えている.

　\(S_W=\Sigma\) と考えて.

t(k)%*%Sigma%*%k

##          LD1
## LD1 1.033366

確かに、\(w^T S_W w=1\)のもとで、.

　　　\({J(w) = \dfrac{w^T S_B w}{w^T S_W w} }=w^T S_B w\) 　　が最大となる (a,b) を求めていることになる.
（wが単位ベクトルでないのは、この理由による）

線形変換されたデータの分散は、元データの分散共分散行列の、線形変換の係数ベクトルによる２次形式に等しい。
つまり、
群内分散を１に保った状態で、群間分散最大の方向 (a,b) を求めていることになる.
実際、
線形変換されたデータの群内分散は１となっている。

YA <- k[1]*dataA[,1]+k[2]*dataA[,2]
YB <- k[1]*dataB[,1]+k[2]*dataB[,2]
(c(var(YA),var(YB)))

## [1] 0.9256546 1.0743454

(c(mean(YA),mean(YB)))

## [1] -0.6100355  3.4698910

mean(c(YA,YB))

## [1] 1.429928

par(mfrow=c(1,2))
hist(YA)
hist(YB)

dist <- data.frame(X=YA-mean(YA),Y=YA-mean(YB))
head(dist)

X	Y
-0.2902920	-4.3702184
0.8804578	-3.1994687
-0.5875778	-4.6675043
3.1892878	-0.8906387
-2.4489766	-6.5289031
-0.8029599	-4.8828864

prob <- data.frame(X=exp(-1/2*(YA-mean(YA))^2),Y=exp(-1/2*(YA-mean(YB))^2))
head(prob)

X	Y
0.9587406	0.0000712
0.6786817	0.0059862
0.8414541	0.0000186
0.0061841	0.6725905
0.0498497	0.0000000
0.7244284	0.0000066

[注意］.
　LD1=k[1] X+ k[2] Y について
LD1はxの線形関数で、ベクトルk方向への(X,Y)の正射影の長さを表すが、各群内でLD1の分散が１となるように係数が定められている。つまりLD1が標準化されているのである。
これは、因子分析でよく行われることである。
LDA関数は LD1が標準化されるように係数を算出している。

群内変動をプールした分散共分散行列は.
\(\begin{split} S_W&=\dfrac{1}{N}\left(\sum_{n \in C_1}(x_n-m_1)(x_n-m_1)'+\sum_{n \in C_2}(x_n-m_2)(x_n-m_2)'\right)\\ &=\dfrac{N_1}{N}\cdot \dfrac{1}{N_1}\sum_{n \in C_1}(x_n-m_1)(x_n-m_1)'+\dfrac{N_2}{N}\cdot \dfrac{1}{N_2}\sum_{n \in C_2}(x_n-m_2)(x_n-m_2)'\\ &=p_1S_1+p_2S_2 \end{split}\).

##          LD1
## X  0.1364599
## Y -0.2488984

t(k)%*%Sigma%*%k

##          LD1
## LD1 1.033366

つまり、
\(w^T S_W w=1\)という条件のもとで \({J(w) = \dfrac{w^T S_B w}{w^T S_W w} }=w^T S_B w\) の最大を考えるので、
ラグランジュの未定乗数法を用いていると考えられる。

　　\(L = w^T S_B w - \lambda (w^T S_W w - 1)\)　を微分して.

　　\(\displaystyle \frac{\partial L}{\partial w} = 2S_B w - 2 \lambda S_W w = 0\)　　\(\therefore S_Bw = \lambda S_W w\).

これから　\((S_W^{-1}S_B) w = \lambda w\)　としているのだと思う。
　　

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

n <- solve(Sigma)%*%(m2-m1)#法線ベクトル
n/norm(n)

##            [,1]
## [1,]  0.3539716
## [2,] -0.6460284

(-n[1]/n[2])

## [1] 0.5479196

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

3.1.1.3 ldaの結果.

## Call:
## lda(Z ~ ., data = d)
## 
## Prior probabilities of groups:
##   A   B 
## 0.5 0.5 
## 
## Group means:
##          X        Y
## A 49.93114 29.82596
## B 69.38844 24.10160
## 
## Coefficients of linear discriminants:
##          LD1
## X  0.1364599
## Y -0.2488984

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

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

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

##        A      B
## 1 0.9999 0.0001
## 2 0.9913 0.0087
## 3 1.0000 0.0000
## 4 0.0091 0.9909
## 5 1.0000 0.0000
## 6 1.0000 0.0000

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

##         LD1
## 1 -2.330255
## 2 -1.159505
## 3 -2.627541
## 4  1.149325
## 5 -4.488940
## 6 -2.842923

var(predict(Z)$x)

##          LD1
## LD1 5.177337

round(mean(predict(Z)$x),1)

## [1] 0

hist(predict(Z)$x)

round(mean(c(YA,YB)),1)

## [1] 1.4

hist(c(YA,YB))

群内分散は１.

var(predict(Z)$x[1:100])

## [1] 0.9256546

var(predict(Z)$x[101:200])

## [1] 1.074345

#プール分散
pv <- (var(predict(Z)$x[1:100])*99+var(predict(Z)$x[101:200])*99)/199

par(mfrow=c(3,1))
hist(predict(Z)$x[1:100],xlim=c(-5,5))
hist(predict(Z)$x[101:200],xlim=c(-5,5))
hist(predict(Z)$x,xlim=c(-5,5))

事後確率. プール分散の値を１として事後確率を計算している.

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

##  [1] 0.382482167 0.270754844 0.335691607 0.002467088 0.019887146 0.289005135
##  [7] 0.335122754 0.226527175 0.295465443 0.337913009

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

##  [1] 2.842211e-05 2.388145e-03 7.416994e-06 2.683248e-01 2.211398e-10
##  [6] 2.651916e-06 7.317950e-06 7.168040e-07 1.694426e-03 8.614996e-04

p1/(p1+p2)

##  [1] 0.999925696 0.991256796 0.999977906 0.009110642 0.999999989 0.999990824
##  [7] 0.999978164 0.999996836 0.994297931 0.997457012

p2/(p1+p2)

##  [1] 7.430412e-05 8.743204e-03 2.209418e-05 9.908894e-01 1.111973e-08
##  [6] 9.175933e-06 2.183615e-05 3.164308e-06 5.702069e-03 2.542988e-03

\(f_q(x)=ax+by-mean(c(YA,YB))\) と考えられる。.

Y <- k[1]*d[,1]+k[2]*d[,2]-mean(c(YA,YB))
round(mean(Y),1)

## [1] 0

head(Y)

## [1] -2.330255 -1.159505 -2.627541  1.149325 -4.488940 -2.842923

これは、LDAのpredict関数で求めた判別関数の値と一致している。

3.1.1.4 事後確率.

線形判別では各群の正規性と等分散性を前提としている。従って、

一次式である判別関数の値は、各群において一次元正規分布に従い.
その分散は各群において等しい.

と考えられる。
つまり、各群の判別関数の値は１次元正規分布と考えられる。
xに対し、

各群の平均からのマハラボビスの距離に基づいて確率が計算され.
条件確率の式から所属確率が算出される.

dnorm() を用いる.

par(mfrow=c(2,1))
hist(predict(Z)$x[1:100],xlim=c(-5,5))
hist(predict(Z)$x[101:200],xlim=c(-5,5))

X <- predict(Z)$x

3.1.1.5 判別関数の値.

２群の正規性と等分散性を仮定しているので、

判別関数の値の集合は 2つの1変量正規分布の混合分布であり.
2つの正規分布の分散は、ともに１である。

２群それぞれの重心は標本平均であるとして、

xが各群に所属する確率を dnorm で計算し.
xのもとでの事後確率を次の式により計算する.
\(\begin{split} P(A|x)&=\dfrac{\pi_Af_A(x)}{\pi_Af_A(x)+\pi_Bf_B(x)} \end{split}\). \(\begin{split} P(B|x)&=\dfrac{\pi_Bf_B(x)}{\pi_Af_A(x)+\pi_Bf_B(x)} \end{split}\).

prob <- data.frame(dnorm(X,mean(X[1:100]),1),dnorm(X,mean(X[101:200]),1))
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

prob <- rename(.data = prob, "A"= "LD1","B"= "LD1.1")
colnames(prob)

## [1] "A" "B"

prob1 <- data.frame(X=prob[1]/(prob[1]+prob[2]),Y=prob[2]/(prob[1]+prob[2]))
head(prob1)

A	B
0.9999257	0.0000743
0.9912568	0.0087432
0.9999779	0.0000221
0.0091106	0.9908894
1.0000000	0.0000000
0.9999908	0.0000092

post <- as.data.frame(predict(Z)$posterior);head(post)

A	B
0.9999257	0.0000743
0.9912568	0.0087432
0.9999779	0.0000221
0.0091106	0.9908894
1.0000000	0.0000000
0.9999908	0.0000092

#f<- post[1]/post[2]
#hist(log(f))
#head(log(f))
head(predict(Z)$posterior)#判別の事後確率

##             A            B
## 1 0.999925696 7.430412e-05
## 2 0.991256796 8.743204e-03
## 3 0.999977906 2.209418e-05
## 4 0.009110642 9.908894e-01
## 5 0.999999989 1.111973e-08
## 6 0.999990824 9.175933e-06

事前確率が異なるときは、さらに平行移動される。
このことを次に見てみよう。

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

A:B=1000:100 とする。

library(MASS) 
mA <- c(50,30)
mB <- c(70,24)
Sigma <- matrix(c(100,35,35,25), 2, 2)
Sigma

##      [,1] [,2]
## [1,]  100   35
## [2,]   35   25

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))

n <- solve(Sigma)%*%(m1-m2);n

##            [,1]
## [1,] -0.5386550
## [2,]  0.9830914

plot(dA[,1:2],xlim=c(20,100),ylim=c(0,50))
par(new=T)
plot(dB[,1:2],xlim=c(20,100),ylim=c(0,50),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] 60.51767 27.19561

m3 <- c(mean(d[,1]),mean(d[,2]))

par(new=T)
plot(x=mean(dA[,1]),y=mean(dA[,2]),xlim=c(20,100),ylim=c(0,50),pch=15,col="blue",xlab="",ylab="")
par(new=T)
plot(x=mean(dB[,1]),y=mean(dB[,2]),xlim=c(20,100),ylim=c(0,50),pch=15,col="blue",xlab="",ylab="")
par(new=T)
plot(m[1],m[2],xlim=c(20,100),ylim=c(0,50),pch=17,col="blue",xlab="",ylab="")
par(new=T)
plot(m3[1],m3[2],xlim=c(20,100),ylim=c(0,50),pch=17,col="red",xlab="",ylab="")

(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 49.88684 29.96341
## B 71.14850 24.42781
## 
## Coefficients of linear discriminants:
##          LD1
## X  0.1392103
## Y -0.2445774

k <- Z $ scaling;k

##          LD1
## X  0.1392103
## Y -0.2445774

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

## [1] -7.25027

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

## [1] 0.5691872

abline(b,a)
abline(b+log(1000/100)/k[2],a,col="blue")

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

## [1] 0.5588945

青い線が事前確率を考慮した場合の判別直線である。

3.1.2.1 母集団比率を考えた線形判別関数.

plot(dA[,1:2],xlim=c(20,100),ylim=c(0,50))
par(new=T)
plot(dB[,1:2],xlim=c(20,100),ylim=c(0,50),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
m3 <- c(mean(d[,1]),mean(d[,2]))

par(new=T)
plot(x=mean(dA[,1]),y=mean(dA[,2]),xlim=c(20,100),ylim=c(0,50),pch=15,col="blue",xlab="",ylab="")
par(new=T)
plot(x=mean(dB[,1]),y=mean(dB[,2]),xlim=c(20,100),ylim=c(0,50),pch=15,col="blue",xlab="",ylab="")
par(new=T)
plot(m3[1],m3[2],xlim=c(20,100),ylim=c(0,50),pch=17,col="red",xlab="",ylab="")

(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 49.88684 29.96341
## B 71.14850 24.42781
## 
## Coefficients of linear discriminants:
##          LD1
## X  0.1392103
## Y -0.2445774

k <- Z $ scaling;k

##          LD1
## X  0.1392103
## Y -0.2445774

b <- m[1]*k[1]/k[2]+m[2]+1/k[2]*log(1000/100);b

## [1] -16.66482

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

## [1] 0.5691872

abline(b,a)

正規母集団で分散共分散行列が等しいが、事前確率が異なる場合、
　　　\(f_q(x)=\log{\dfrac{\pi_1}{\pi_2}}-\dfrac{1}{2}D_1^2+\dfrac{1}{2}D_2^2\)　　が線形判別関数となる。

head(predict(Z)$x)

##          LD1
## 1 -1.5942848
## 2 -1.2983494
## 3 -0.1544567
## 4 -1.6176139
## 5 -0.1793112
## 6 -1.3379446

log(100/1000)

## [1] -2.302585

head(log(100/1000)+predict(Z)$x)

##         LD1
## 1 -3.896870
## 2 -3.600934
## 3 -2.457042
## 4 -3.920199
## 5 -2.481896
## 6 -3.640530

Y <- k[1]*d[,1]+k[2]*d[,2]-mean(c(YA,YB))
round(mean(Y),1)

## [1] -1.4

t(k)%*%Sigma%*%k

##          LD1
## LD1 1.050065

YA <- k[1]*dataA[,1]+k[2]*dataA[,2]
YB <- k[1]*dataB[,1]+k[2]*dataB[,2]
(c(var(YA),var(YB)))

## [1] 1.0307949 0.6892513

Y <- k[1]*d[,1]+k[2]*d[,2]-mean(c(YA,YB))+log(100/1000)
round(mean(Y),1)

## [1] -2.3

head(Y)

## [1] -3.896870 -3.600934 -2.457042 -3.920199 -2.481896 -3.640530

事前確率による修正項 log(100/1000) を加えた式が判別関数として使われていることがわかる.

3.1.3 事後確率.

次の式によって事後確率を算出し、LDAの出力と比べる。

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

head(predict(Z)$x)

##          LD1
## 1 -1.5942848
## 2 -1.2983494
## 3 -0.1544567
## 4 -1.6176139
## 5 -0.1793112
## 6 -1.3379446

round(mean(predict(Z)$x),1)

## [1] 0

head(round(predict(Z)$posterior,5))

##         A     B
## 1 1.00000 0e+00
## 2 1.00000 0e+00
## 3 0.99997 3e-05
## 4 1.00000 0e+00
## 5 0.99998 2e-05
## 6 1.00000 0e+00

hist(predict(Z)$x)

par(mfrow=c(2,1))
hist(predict(Z)$x[1:1000],xlim=c(-5,5))
hist(predict(Z)$x[1001:1100],xlim=c(-5,5))

X <- predict(Z)$x
prob1 <- dnorm(X,mean(X[1:1000]),1)
prob2 <- dnorm(X,mean(X[1001:1100]),1)
prob <- data.frame(dnorm(X,mean(X[1:1000]),1),dnorm(X,mean(X[1001:1100]),1))
prob <- rename(.data = prob, "A"= "LD1","B"= "LD1.1")
colnames(prob)

## [1] "A" "B"

prob1 <- data.frame(X=1000/1100*prob[1]/(1000/1100*prob[1]+100/1100*prob[2]),Y=100/1100*prob[2]/(1000/1100*prob[1]+100/1100*prob[2]))
head(prob1)

A	B
0.9999999	1.00e-07
0.9999998	2.00e-07
0.9999746	2.54e-05
1.0000000	0.00e+00
0.9999772	2.28e-05
0.9999998	2.00e-07

post <- as.data.frame(predict(Z)$posterior);head(post)

A	B
0.9999999	1.00e-07
0.9999998	2.00e-07
0.9999746	2.54e-05
1.0000000	0.00e+00
0.9999772	2.28e-05
0.9999998	2.00e-07

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

##           A            B
## 1 0.9999999 5.095415e-08
## 2 0.9999998 1.826387e-07
## 3 0.9999746 2.538520e-05
## 4 1.0000000 4.607594e-08
## 5 0.9999772 2.280440e-05
## 6 0.9999998 1.539622e-07

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

A:B=1000:500 とする。

library(MASS) 
mA <- c(50,30)
mB <- c(70,24)
Sigma <- matrix(c(100,35,35,25), 2, 2)
Sigma

##      [,1] [,2]
## [1,]  100   35
## [2,]   35   25

dataA <- mvrnorm(1000, mA, Sigma)
dataB <- mvrnorm(500, 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",500))

n <- solve(Sigma)%*%(m1-m2);n

##            [,1]
## [1,] -0.5688528
## [2,]  1.0178177

plot(dA[,1:2],xlim=c(20,100),ylim=c(0,50))
par(new=T)
plot(dB[,1:2],xlim=c(20,100),ylim=c(0,50),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] 60.03071 26.93830

m3 <- c(mean(d[,1]),mean(d[,2]))

par(new=T)
plot(x=mean(dA[,1]),y=mean(dA[,2]),xlim=c(20,100),ylim=c(0,50),pch=15,col="blue",xlab="",ylab="")
par(new=T)
plot(x=mean(dB[,1]),y=mean(dB[,2]),xlim=c(20,100),ylim=c(0,50),pch=15,col="blue",xlab="",ylab="")
par(new=T)
plot(m[1],m[2],xlim=c(20,100),ylim=c(0,50),pch=17,col="blue",xlab="",ylab="")
par(new=T)
plot(m3[1],m3[2],xlim=c(20,100),ylim=c(0,50),pch=17,col="red",xlab="",ylab="")

(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 49.75103 29.84266
## B 70.31040 24.03393
## 
## Coefficients of linear discriminants:
##          LD1
## X  0.1384401
## Y -0.2428012

k <- Z $ scaling;k

##          LD1
## X  0.1384401
## Y -0.2428012

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

## [1] -7.289939

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

## [1] 0.5701788

abline(b,a)
abline(b+log(1000/500)/k[2],a,col="blue")

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

## [1] 0.5515701

注意.

法線ベクトルの理論ベクトルとLDAのscalingで求めたベクトルは符号が逆になっていることに注意。このことは、事前確率による切片の修正をする場合の符号に影響する。

事前確率が大きく異なれば、判別直線は移動する。

判別分析のまとめ