#はじめに
本記事は『The Elements of Statistical Learning』のゼミ用記事です。

#発表範囲
4.Linear Methods for Classification(P101~)
-4.1. Introduction
-4.2. Linear Regression of an Indicator Matrix
-4.3. Linear Discriminant Analysis

#参考文献
分類の為の線形手法(@kenmatsu4 )

#4.1 Introduction

これまでの章で我々は主に回帰問題を扱ってきたが、この章では分類問題について考える。
2章で分類問題の解き方について簡単に触れた（k-means法と最小二乗法）

ここでクラス分けのアプローチとしては
1.判別関数　$\delta_k(x)$　が最大となるクラス$k$
2.事後分布 ${\rm Pr}(G=k|X=x)$　が最大となるクラス$k$

という2種類のパターンが考えられる。

##判別関数
$k$番目の指示応答変数に対する線形モデル

\delta_k(x) = \hat\beta_{k0} + \hat\beta^T_kx

で表される。

この時クラス$k$と$l$の決定境界は

\{ x:(\hat\beta_{k0} - \hat\beta_{l0}) + (\hat\beta_{k} - \hat\beta_{l})^Tx = 0 \}

でありこれは超平面となる。（一般に超平面は原点を通る場合の事をさすが、この章ではアフィン集合であっても超平面と言う事にする）

またここで任意のクラスペアで決定境界が超平面である事が明らかなので、入力空間は区分的に超平面によって分割する事ができる。

##事後分布
ここで事後分布として以下のようなものを考える。

$${\rm Pr}(G=1|X=x) = \frac{exp(\beta_0+\beta^Tx)}{1+exp(\beta_0+\beta^Tx)}$$

$${\rm Pr}(G=2|X=x) = \frac{1}{1+exp(\beta_0+\beta^Tx)}$$

これをlogit変換:$log[p/(1-p)]$ をすると

$${\rm log}\frac{{\rm Pr}(G=1|X=x)}{{\rm Pr}(G=2|X=x)} = \beta_0+\beta^Tx $$

というlog-odds関数を得る事ができる。

この時決定境界は

 \{ x|\beta_0+\beta^Tx=0 \}

で表す事ができる。

##拡張された入力空間を扱う

本章は線形な決定境界に関して議論するが、一般化して考える際には注意が必要である。
入力空間が2次元の場合の入力変数のセットが$X_1,X_2$の時は2次元空間で決定境界が線形になるが、$X_1, X_2, X_1^2, X_2^2, X_1X_2$の5次元の線形関数にフィットさせたものを2次元に映すとそれは曲線を表す。

#4.2 Linear Regression of an Indicator Matrix

前章で多変量のアウトプットを扱った問題を取り上げた。
今回はそれを利用して分類問題を考えてみる。

一般にカテゴリが$K$個あった時、出力$Y\in \mathbb{R}^{N\times K }$、入力$X \in \mathbb{R}^{N\times (p+1)}$の行列である

また$Y$の各行は正解ラベルのみが１になっていて残りはすべて０である。

この時線形回帰手法の推定値は

\hat{Y} = X(X^TX)^{-1}X^TY
\\
\hat{B} = (X^TX)^{-1}X^TY

と求める事ができる。

ここで新たなデータ点$x$が与えられたとしたら推定値$\hat B$を用いて予測値を出し、値が高かったラベルに分類する。

\hat{f}(x) = (1 ,x^T)\hat{B}　\\
\hat{G}(x) = {\rm argmax}_{k\in \mathcal{G}}\hat{f}_k(x)

この手法の合理的解釈として$\hat f(x)$が確率${\rm Pr}(G=k|X=x)$として考えられればとして考えられればある程度妥当な判断と言えるだろう。

ただし、$\sum_{k\in\mathcal{G}}\hat{f}_k(x) = 1$を満たさない事や$\hat{f}_k(x)$がマイナスの値を取ったり1より大きくなってしまう事はよくある（主に訓練データ外で）

実際には基底展開してから線形回帰を行う事で確率分布に一致推定量を得る事ができる。（詳細は5章）

ここでより単純な見方をして、上の線形回帰モデルを最小二乗法の線形モデルとして考えると

\underset{B}{{\rm min}}\sum_{i=1}^N||y_i-[(1,x^T_i)B]^T||^2 \\

\hat{G}(x) = \underset{k}{\rm argmin}||\hat{f}(x)-t_k||^2

と解釈できる。

ここで$t_k$は$K\times K$の単位行列の$k$列目を表す。

線形回帰モデルは$K\geq3$の時は上手く分離できない事が多くある。
以下は$K=3$の時の線形回帰モデルでの分割である（左図）
右図はLDAによる分割（後述）

このように線形回帰モデルでの分類ではクラスが3個以上になった時に上手く分類できなくなってしまう事が多くなる。
下の左図に示したように、入力を$X_1$と$X_2$の2入力の空間にした時には、青色のクラスが他のクラスの判別関数にかぶせられてしまい青色のクラスがマスクされてしまう。
ただし、入力空間を拡張して二乗項まで考えた時は上手く分類できている事が右図より分かる。

一般には$K$個のクラス分けの場合は$K-1$次まで考慮したほうが上手く分類できる可能性が高いと考えられる。

また、入力の$p$の次元数よりもクラス$K$の数のほうが多い場合は上手く分類することが難しい事が知られている。

線形モデルの実装（Python)

import numpy as np
import matplotlib.pyplot as plt

x1 = np.random.normal(10,3,100)
y1 = np.random.normal(10,3,100)
x2 = np.random.normal(0,3,100)
y2 = np.random.normal(0,3,100)
x3 = np.random.normal(-10,3,100)
y3 = np.random.normal(-10,3,100)

X1 = np.c_[x1,y1]
X2 = np.c_[x2,y2]
X3 = np.c_[x3,y3]

data = np.vstack([X1,X2,X3])
data = np.c_[np.ones(data.shape[0]),data]

lab1 = np.array([0,0,1]).reshape(1,-1).repeat(100,axis=0)
lab2 = np.array([0,1,0]).reshape(1,-1).repeat(100,axis=0)
lab3 = np.array([1,0,0]).reshape(1,-1).repeat(100,axis=0)
label = np.r_[lab1,lab2,lab3]

B = np.dot(np.dot(np.linalg.inv(np.dot(data.T,data)),data.T),label)
grid = np.meshgrid(np.linspace(-20,  20, 401), np.linspace(-20, 20, 401))
vel = np.c_[np.ones(401*401),grid[0].flatten(),grid[1].flatten()]
z = np.argmax(np.dot(vel,B),axis=1).reshape(401,401)

plt.figure(figsize=(10,10))
plt.pcolor(grid[0], grid[1], z, alpha=0.4)
plt.scatter(x1,y1,color='r')
plt.scatter(x2,y2,color='g')
plt.scatter(x3,y3,color='b')
plt.xlim(-20,20)
plt.ylim(-20,20)

x = np.linspace(-20,20,401)
y = np.ones(401)
field = np.c_[np.ones(401), np.linspace(-20,20,401), y]
ans = np.dot(field,B)

plt.plot(x,ans[:,0],color='b')
plt.plot(x,ans[:,1],color='g')
plt.plot(x,ans[:,2],color='r')

線形モデルの実装（拡張空間）（Python)

x1 = np.random.normal(10,3,100)
y1 = np.random.normal(10,3,100)
x2 = np.random.normal(0,3,100)
y2 = np.random.normal(0,3,100)
x3 = np.random.normal(-10,3,100)
y3 = np.random.normal(-10,3,100)

X1 = np.c_[x1,y1,x1**2,y1**2,x1*y1]
X2 = np.c_[x2,y2,x2**2,y2**2,x2*y2]
X3 = np.c_[x3,y3,x3**2,y3**2,x3*y3]

data = np.vstack([X1,X2,X3])
data = np.c_[np.ones(data.shape[0]),data]

lab1 = np.array([0,0,1]).reshape(1,-1).repeat(100,axis=0)
lab2 = np.array([0,1,0]).reshape(1,-1).repeat(100,axis=0)
lab3 = np.array([1,0,0]).reshape(1,-1).repeat(100,axis=0)
label = np.r_[lab1,lab2,lab3]

B = np.dot(np.dot(np.linalg.inv(np.dot(data.T,data)),data.T),label)

grid = np.meshgrid(np.linspace(-20,  20, 401), np.linspace(-20, 20, 401))

vel = np.c_[np.ones(401*401),grid[0].flatten(),grid[1].flatten(),grid[0].flatten()**2,grid[1].flatten()**2,(grid[0]*grid[1]).flatten()]

z = np.argmax(np.dot(vel,B),axis=1).reshape(401,401)

plt.figure(figsize=(10,10))
plt.pcolor(grid[0], grid[1], z, alpha=1)
plt.scatter(x1,y1,color='r')
plt.scatter(x2,y2,color='g')
plt.scatter(x3,y3,color='b')
plt.xlim(-20,20)
plt.ylim(-20,20)


x = np.linspace(-20,20,401)
field = np.c_[np.ones(401), np.linspace(-20,20,401), np.zeros(401),x**2, np.zeros(401), np.zeros(401)]
ans = np.dot(field,B)

plt.plot(x,ans[:,0],color='b')
plt.plot(x,ans[:,1],color='g')
plt.plot(x,ans[:,2],color='r')

#4.3 Linear Discriminant Analysis

分類問題においては事後分布${\rm Pr}(G|X)$を求める事が大きな目標になる。（詳細：Section 2.4）

以降で出てくる式を列挙しておく。

f_k(x):データXの条件付き確率 \rightarrow {\rm Pr}(X=x|G=k)　\\ 
\pi_k:カテゴリの事前分布 \rightarrow {\rm Pr}(G=k)

ここで、ベイズの定理より、

{\rm Pr}(G=k|X=x) = \frac{f_k(x)\pi_k}{\sum^{K}_{l=1}f_l(x)\pi_l}

が成り立つ。

各クラスのデータの分布については様々なケースが考えられるが、今回は多変量ガウス分布であるケースを想定する。
(その他の分布についてはSection 6で扱う。)

f_k(x) = \frac{1}{(2\pi)^{p/2} |\Sigma_k|^{1/2}}{\rm exp}(-\frac{1}{2}(x-\mu_k)^T\Sigma^{-1}_k(x-\mu_k))

Linear discriminant analysis(LDA)は分散共分散行列$\Sigma_k$が全て同じ$\Sigma$であるという仮定を置いて考える。

\Sigma_k = \Sigma \quad \forall k

この時クラス$k$と$l$の事後分布の$\rm log$比を見てみると、

\begin{align}
{\rm log}\frac{{\rm Pr}(G=k|X=x)}{{\rm Pr}(G=l|X=x)} &= {\rm log}\frac{f_k(x)}{f_l(x)}+{\rm log}\frac{\pi_k}{\pi_l}
\\
&={\rm log}\frac{\pi_k}{\pi_l}-\frac{1}{2}(\mu_k+\mu_l)^T\Sigma^{-1}(\mu_k-\mu_l) + x^T\Sigma^{-1}(\mu_k-\mu_l)

\end{align}

を得る事ができ、　$x$の線形関数になっている事が分かる。

${\rm Pr}(G=k|X=x)={\rm Pr}(G=l|X=x)$の時がクラス$k$と$l$の決定境界となる。

この手法での決定境界は垂直二等分線にはなっていない事が特徴になっている。

ここで判別関数は以下のようになる。

\delta_k(x) = x^T\Sigma^{-1}\mu_k - \frac{1}{2}\mu_k^T\Sigma^{-1}\mu_k + {\rm log}\pi_k

この判別関数は以下のルールで分類する事で事後確率での議論と等価になる。

G(x) = {\rm argmax}_k\delta_k(x)

実問題として、分布の平均値や事前分布は分からないので以下の方法で推定する。

\hat\pi_k = N_k/N 
\\
\hat\mu_k = \sum_{g_i=k}x_i/N_k 
\\
\hat\Sigma = \sum^K_{k=1}\sum_{g_i=k}(x_i-\hat\mu_k)(x_i-\mu_k)^T/(N-K)

2クラス分類の時、LDAを用いると以下の不等式を満たせばクラス２に分類され、それ以外の時は１に分類される。

x^T\hat\Sigma^{-1}(\hat\mu_2-\hat\mu_1) > \frac{1}{2}\hat\mu_2^T\Sigma^{-1}\hat\mu_2 - \frac{1}{2}\hat\mu_1^T\Sigma^{-1}\hat\mu_1 + {\rm log}(N_1/N) - {\rm log}(N_2/N)

##異なる分散共分散行列を持つ場合
ここでLDAでは$\Sigma_k$が全て共通であるという仮定を置いていたが、各クラスによって異なる分散共分散行列を持つ場合、分類関数は

\delta_k(x) = -\frac{1}{2}{\rm log}|\Sigma_k| - \frac{1}{2}(x-\mu_k)^T\Sigma^{-1}_k(x-\mu_k) + {\rm log}\pi_k

となる。クラス$k$と$l$の決定境界は

\{x:\delta_k(x) = \delta_l(x) \}

で表される。

QDAは各分散共分散行列を推定しなければいけない

下図で示す通り、LDAとQDAは大きく差はないが、一般的にQDAの方が誤差が少ない。

左図は入力空間を5次元に拡張した場合のLDA、右図は2次元のQDAである。

ただし、次元が大きくなるにつれてQDAのパラメータ数はとても大きくなる。

	パラメータ数
LDA	$(K-1)\times (p+1)$
QDA	$(K-1)\times ( p(p+3)/2+1 ) $

その為入力空間の次元が極端に大きい場合などはLDAで代用する事が好まれる事もある。

STATLOGプロジェクト(Michie et al., 1994)にて22のデータセットに対して様々な分類手法を用いてその性能の良さを調べた所
LDAは22個中7個のデータセットでトップ3に入っており、QDAは22個中4つがトップ3に入っています。
10個のデータセットにはLDAかQDAのいずれかがトップ3に入っているという結果が示されました。

LDAとQDAが非常にデータに対して有効である原因としてはデータは多くがガウス分布によって生成されていると考えられる。