フィッシャーの線形判別分析法

# はじめに
　学校課題のついでに，線形判別分析（Linear Discriminant Analysis, LDA）の有名なアルゴリズムであるFisherの線形判別について書いてみました．分かりにくい部分もあると思いますが，ご容赦ください．

# 線形判別分析
　線形判別分析とは2つのクラスを"**最もよく判別できる**"直線を求める手法です．
　データが直線のどちら側にあるかを見ることで，どちらのクラスに属するか判別することができます．この直線のことを**決定境界**と言います．わかりやすいように各クラスに属するデータは2次元のデータであるとすると，決定境界は下図のようなイメージです．

<img src="https://qiita-image-store.s3.amazonaws.com/0/187083/bd247f4b-1f8d-1e4d-7da3-f5561b13c247.png" width=60%>

　「直線のどちら側にあるかを見る」と言いましたが，実際にデータを識別するときは適当な線形変換によってデータを直線<i>z</i>上に写像すると2つのクラスを判別しやすいと思います．データを1次元に圧縮すれば，閾値となるスカラー値を適当に設定してどちらのクラスに属するのか判別できますね．フィッシャーの線形判別分析は，この射影先の直線<i>z</i>を見つけるための手法の1つです．
　
<img src="https://qiita-image-store.s3.amazonaws.com/0/187083/ac57e353-d097-8196-70a0-c0ab8cca2d6c.png" width=60%>

　さて，この線形変換を行列式で表すと，

```math
y_n = w^T x_n
```

となります．<i>x<sub>n</sub></i>はD次元のデータ，<i>w</i>が線形変換ベクトル（<i>x<sub>n</sub></i>と同じD次元のベクトル），<i>y<sub>n</sub></i>は線形変換された1次元データ（スカラー）です．

　ちなみに線形判別分析は，判別対象が2クラスではなく多クラスになっても対応可能です．
　また，各クラスのデータが2次元ではなく3次元以上でも適応できます．この場合の決定境界は直線ではなく平面や超平面になります．

# フィッシャーの線形判別分析の導出
　実際にデータを見ながらフィッシャーの線形判別の導出をやってみます．フィッシャーの線形判別は，<font color=#ff0000>射影した各クラスの分離度が最も大きくなるような線形変換ベクトル<i>w</i>（データを射影する直線<i>z</i>の傾きになる）を見つける</font>という手法になります．

　まず，以下のようなデータ<i>C<sub>1</sub></i>と<i>C<sub>2</sub></i>を考えます．

<img src="https://qiita-image-store.s3.amazonaws.com/0/187083/23a5b9bd-08e1-4724-efd5-7291426aeca3.png" width=60%>

　クラス<i>C<sub>k</sub></i>のデータ数を<i>N<sub>k</sub></i>とすると，各クラスの平均は，

```math
m_1 = \frac{1}{N_1}\sum_{n \in C_1}x_n
```

```math
m_2 = \frac{1}{N_2}\sum_{n \in C_2}x_n
```

となります．この<i>m<sub>1</sub></i>と<i>m<sub>2</sub></i>を各クラスの代表点として，この2点が線形変換ベクトル<i>w</i>によって射影されたときに最も離れるようにします．つまり，

```math
w^T m_1 - w^T m_2 = w^T (m_1 - m_2)
```

が最大になるように<i>w</i>を決定します．ただし，必要なのは最大値を与えるベクトル<i>w</i>の傾きなのでベクトル<i>w</i>自体の大きさを大きくしても意味がありません．したがって，<i>w</i>は単位長であるという制約

```math
|w|=1
```

を与えて，ラグランジュの未定乗数法で最大値を求めると

```math
w \propto (m_1 - m_2)
```

となります．これは，<i>w</i>が各クラスの平均値を結ぶ線分と平行になるということを意味しています．実際にこのベクトル<i>w</i>から決定境界の直線を引いてみます．決定境界の直線はwに直交し，各クラスの平均を結ぶ線分の中点を通る直線になっているので以下の図のようになります．

<img src="https://qiita-image-store.s3.amazonaws.com/0/187083/7e3ba824-0f22-b98c-4550-3750c1ac3de3.png" width=60%>

<i>C<sub>1</sub></i>の青い点が決定境界の直線の上にもあります．これだけでは上手く判別することができませんでした．

***

　実は，各クラスの平均値だけを見ていると上手く判別することができません．例えば，<i>C<sub>1</sub></i>の青い点が<i>C<sub>2</sub></i>の赤い点が密集している位置に混ざっていても，一方で<i>C<sub>2</sub></i>から遠く離れた位置にたくさんの青い点が存在すれば，平均的には2つのクラスは離れているとみなされてしまうからです．つまり，各クラスの点のばらつきが大きいと，平均的な位置が離れていても適切な決定境界を引くことができないということです．

　ここでベクトル<i>w</i>の決定のためにもう1つ条件を加えます．それは，同じクラス内のデータは射影した後でなるべく近くに密集してほしいというものです．つまり，射影後の各クラス内の分散を小さくするということです．射影後の各クラスの分散を数式で表すと，

```math
s_1^2 = \sum_{n \in C_1} (w^T x_n - w^T m_1)^2
```

```math
s_2^2 = \sum_{n \in C_2} (w^T x_n - w^T m_2)^2
```

となります．したがって，各クラス内の分散を小さくするためには，

```math
s^2 = s_1^2 + s_2^2
```

を最小化すれば良いということです．

***

 適切な<i>w</i>を見つけるためには，

1. 線形変換した2つのクラスがなるべく離れるように，各クラス間の平均の差がなるべく大きくなるようにする．
2. 線形変換した各クラス内のデータがなるべく密集するように，各クラス内の分散がなるべく小さくなるようにする．

の2つの条件が必要なことがわかりました．この2つの条件は，分数を用いて1つの条件として扱うことができます．つまり，

```math
J(w) = \frac{(クラス間の平均に関する最大化問題)}{(クラス内の分散に関する最小化問題)}
```

という式にして，評価関数<i>J(w)</i>を最大化すれば良いということです．

　ただし，分子と分母を同様に扱うために，分子も平均ではなく分散の形で表します.
つまり分子を，

```math
w^T (m_1 -m_2)
```

から

```math
w^T (m_1 - m_2) (m_1 - m_2)^T w
```

という形式にします．このとき，

```math
S_B = (m_1 - m_2) (m_1 - m_2)^T
```

を**クラス間共分散行列**と呼びます．

したがって，分子は

```math
w^T S_B w
```

と書くことができます．

　また，分母は，

```math
s^2 = s_1^2 + s_2^2 = \sum_{n \in C_1} (w^T x_n - w^T m_1)^2 + \sum_{n \in C_2} (w^T x_n - w^T m_2)^2
```

でしたが，

```math
S_W = \sum_{n \in C_1} (x_n - m_1)^2 + \sum_{n \in C_2} (x_n - m_2)^2
```

とおくと，

```math
s^2 = w^T S_W w
```

と書くことができます．このとき<i>S<sub>W</sub></i>を**総クラス内共分散行列**と呼びます．

　したがって，評価関数<i>J(w)</i>は，

```math
J(w) = \frac{w^T S_B w}{w^T S_W w}
```

と表現することができて，これを最大化する<i>w</w>を求めます．ラグランジュの未定乗数法で解くと最終的に，

```math
w \propto {S_W}^{-1} S_B w \propto {S_W}^{-1} (m_1 - m_2)
```

　となります．

　結論は，D次元のデータから各クラス内の平均と総クラス内共分散行列の逆行列を求めれば線形変換ベクトル<i>w</i>が求まり，決定境界が定まるということです．

<img src="https://qiita-image-store.s3.amazonaws.com/0/187083/4e6f0033-d06e-7806-d178-55c946f73b79.png" width=60%>

実際にやってみると，良い感じに決定境界の直線を引くことができました．

***

 適当にデータを生成して，フィッシャーの線形判別で決定境界を引くプログラムを書いてみる．

```fisher.py
#coding: utf-8

import numpy as np
import matplotlib.pyplot as plt

DIM = 2 #データの次元数

#直線の式
def f(x, a, b):
	return a*x+b

#図示
def plot(cls1, cls2, line=None):
	x,y = cls1.T
	plt.plot(x, y, 'bo', ms=3, label='class1')
	x, y = cls2.T
	plt.plot(x, y, 'ro', ms=3, label='class2')

	if not (line is None):
		plt.plot(line[0], line[1], 'k-', ms=5)

	plt.xlim(-10,10)
	plt.ylim(-10,10)

	plt.show()

def fisher(cls1, cls2):
	#リストからnp.arrayに変換（行列の転置や逆行列を扱うため）
	cls1 = np.array(cls1)
	cls2 = np.array(cls2)

	#各クラスの平均値
	mean1 = np.mean(cls1, axis=0)
	mean2 = np.mean(cls2, axis=0)

	#総クラス内共分散行列
	sw = np.zeros((DIM,DIM))
	for xn in cls1:
		xn = xn.reshape(DIM,1)
		mean1 = mean1.reshape(DIM,1)
		sw += np.dot((xn-mean1),(xn-mean1).T)
	for xn in cls2:
		xn = xn.reshape(DIM,1)
		mean2 = mean2.reshape(DIM,1)
		sw += np.dot((xn-mean2),(xn-mean2).T)

	#総クラス内共分散行列の逆行列
	sw_inv = np.linalg.inv(sw)

	#wを求める
	w = np.dot(sw_inv,(mean1-mean2))

	#決定境界直線を図示する
	mean = (mean1 + mean2)/2 #平均値の中点
	a = -w[0]/w[1] #wと直交する
	b = a*mean[0]+mean[1]
	x = np.linspace(-8, 8, 1000)
	y = f(x,a,b)

	plot(cls1, cls2, (x,y))
	
	return w

if __name__ == '__main__':
	#テスト用2次元データ
	cov = [[3,1],[1,3]] #共分散
	cls1 = np.random.multivariate_normal([-5,-5], cov, 50)
	cls2 = np.random.multivariate_normal([5,5], cov, 50)

	#fisherで線形変換ベクトルwを求める
	w = fisher(cls1, cls2)
	print('w =', w, '\n')
```

実行結果

<img src="https://qiita-image-store.s3.amazonaws.com/0/187083/aa8ef346-d971-0c01-3d43-4399a542c049.png" width=60%>

良い感じの決定境界の直線が引けてますね．

# おわりに
　Fisherの線形判別分析をやってみました．他の手法との識別精度の比較をしてみたいですね．

# 参考
* [フィッシャーの線形判別](http://s0sem0y.hatenablog.com/entry/2016/08/07/054847)
* [フィッシャーの線形判別](http://aidiary.hatenablog.com/entry/20100425/1272158587)
* [線形識別器 - Fisherの線形判別関数，単純パーセプトロン，IRLS](https://qiita.com/39Cu4c/items/be8ace83dc958b124921)