PRML-4.1 判別函式
1.概念
判別式是一個使用輸入向量\(x\)並把它分配給\(K\)種分類的其中一種\(C_k\)的函式。本章中，我們把我們的討論侷限於線性判別式（linear discriminants），即那些決策面是超平面的判別函式。為了簡化討論，我們首先考慮二分類的情況，再推廣到\(K > 2\)的情形。
2 二分類
線性判別式的最簡單形式是取輸入向量的線性函式，即
$ y(x) = w^Tx + w_0 \tag{4.4} $
如果\(y(x) \geq 0\)則把輸入向量\(x\)分到類\(C_1\)中，否則分到\(C_2\)中。因此，對應的決策邊界由\(y(x) = 0\)確定，它對應著\(D\)
原點到決策面的標準距離由
$ \frac{w^Tx}{\Vert w \Vert} = -\frac{w_0}{\Vert w \Vert} \tag{4.5} $
解釋一下
設點\(x\)在決策面上,\(w\)是法向量
\(\frac{w^Tx}{\Vert w \Vert} 代表向量x在w上的投影，=\frac{||w||||x||cos\theta}{\Vert w \Vert}=||x||cos\theta\)
因為x在決策面上,\(w^Tx + w_0=0 ,w^Tx = - w_0\)
3 多分類
現在把線性判別式推廣到\(K > 2\)
考慮一個包含\(K-1\)個把屬於類別\(C_k\)的點與其它的區分開的分類器(本類和其他類)。這被稱為一對其他（one-versus-the-rest）分類器。圖4.2的左手邊展示了一個涉及三個類別的例子，這種方法會導致空間中有一塊類別含糊不清的區域。
另一種方法是引入\(K(K-1)/2\)個二元判別函式，對每一對類別都設定一個。這被稱為一對一（one-versus-one）分類器。然後根據這些判別函式的大多數投票來確定每個點的類別。然而這也會引起類別含糊不清的區域的問題，圖4.2的右手邊展示。
所以，引入包含\(K\)個形式為
\(y_k(x) = w_k^Tx + w_{k0} \tag{4.9}\)
的線性函式的單個\(K\)類判別式，並把點\(x\)分入對於所有\(j \neq k\)都有\(y_k(x) > y_j(x)\)的類\(C_k\)中，就可以避免這些問題。類\(C_k, C_j\)間的決策邊界由\(y_k(x) = y_j(x)\)給出，它對應形式為
$ (w_k - w_j)^Tx + (w_{k0} - w_{j0}) = 0 \tag{4.10} $
的\((D-1)\)維的超平面。這與4.1.1節討論的二分類情況下的決策邊界具有相同的形式，因此也有類似的幾何性質。
這樣的判別式得到的決策區域總是單連通且凸的。為了證明這個，讓我們考慮圖4.3中展示的，兩個都位於決策區域\(R_k\)中的點\(x_A,x_B\)。
任何一個在連線\(x_A,x_B\)線段上的點\(\hat{x}\)可以表示為
$ \hat{x}=\lambda x_A + (1-\lambda)x_B \tag{4.11} $
其中\(0 \leq \lambda \leq 1\)。根據判別函式的線性性，可以得到
$ y_k(\hat{x}) = \lambda y_k(x_A) + (1-\lambda)y_k(x_B) \tag{4.12} $
因為\(x_A, x_B\)都在\(R_k\)中，所以對於所有\(j \neq k\)都滿足\(y_k(x_A) > y_j(x_A), y_k(x_B) > y_j(x_B)\)，所以\(y_k(\hat{x}) > y_j(\hat{x})\)，所以\(\bar{x}\)在\(R_k\)中。因此\(R_k\)是單連通且凸的(我理解是如果兩個點都屬於第k類，它們之間的點也都是屬於k類別)。
注意，對於二分類的情形，我們既可以使用這裡討論的基於兩個判別函式\(y_1(x), y_2(x)\)方法，也可以使用4.1.1節給出基於單一的判別函式\(y(x)\)的更簡單的但等價的方法。
現在，我們開始探討三種線性判別函式的引數學習方法，即基於最小二乘、Fisher線性判別式，以及感知器演算法。
4.判別函式1-最小平方法法
每個類別\(C_k\)通過它自己的線性模型
$ y_k(x) = w_k^Tx + w_{k0} \tag{4.13} $
可以很容易地把這些量聚集在一起表示，即
$ y(x) = \widetilde{W}^T\tilde{x} \tag{4.14} $
其中
\[\widetilde{W} 是第 k^{th} 列構成 D+1 維向量 \tilde{w}k = (w_{k0}, w_k^T)^T 的矩陣 \]\(\tilde{x}\)
\[\tilde{x} 對應的是通過虛輸入x_0=1擴張後的輸入向量(1,x^T)^T \]訓練資料集
\[{x_n,t_n} n = 1,...,N \]\(矩陣T\),y_test
\[定義一個第 n^{th} 行是向量t_n^T 矩陣T \]\(矩陣X\),x_test
\[行是\tilde{x}_n^T的矩陣\tilde{X} \]平方和誤差函式就可以寫成
\[E_{D}(\tilde{\boldsymbol{W}})=\frac{1}{2} \operatorname{Tr}\left\{(\tilde{\boldsymbol{X}} \tilde{\boldsymbol{W}}-\boldsymbol{T})^{T}(\tilde{\boldsymbol{X}} \tilde{\boldsymbol{W}}-\boldsymbol{T})\right\}\tag{4.15} \]證明 4.15
先詳細展開一些標記
\(\tilde x =\begin{pmatrix}
1 \\
x_1\\
.\\
x_D\\
\end{pmatrix}_{D+1 \times 1},\tilde X =\begin{pmatrix}
\tilde x_1^T\\
.\\
\tilde x_N^T\\
\end{pmatrix}_{N \times D+1},T =\begin{pmatrix}
t_1^T \\
t_2^T\\
.\\
t_N^T\\
\end{pmatrix}_{N\times K},\tilde W = \begin{pmatrix}
w_{10} & ... & w_{K0} \\
... & ... & ... \\
w_{1D} & ... & w_{KD} \\
\end{pmatrix}_{D+1 \times K},每一列都是\begin{pmatrix}
w_{k0} \\
... \\
w_{kD} \\
\end{pmatrix} = \tilde w_k,\tilde W^T = \begin{pmatrix}
w_{10} & ... & w_{1D} \\
... & ... & ... \\
w_{K0} & ... & w_{KD} \\
\end{pmatrix}_{K \times D+1}\)
\(y(x) = \widetilde{W}^T\tilde{x}= \begin{pmatrix}
w_{10} & ... & w_{1D} \\
... & ... & ... \\
w_{K0} & ... & w_{KD} \\
\end{pmatrix}_{K \times D+1}\begin{pmatrix}
1 \\
x_1\\
.
x_D\\
\end{pmatrix}_{ D+1 \times 1}\)
\(y(x)=\begin{pmatrix}
y_1(x) \\
... \\
y_K(x) \\
\end{pmatrix}_{K\times 1}=\begin{pmatrix}
p_1(x) & 表示資料x歸類於類別1的概率\\
... \\
p_K(x) & 表示資料x歸類於類別K的概率\\
\end{pmatrix}_{K\times 1}\)
然後開始計算
\(\tilde X \tilde W = \begin{pmatrix}
\tilde x_1^T\\
.\\
\tilde x_N^T\\
\end{pmatrix}\begin{pmatrix}
w_{10} & ... & w_{K0} \\
... & ... & ... \\
w_{1D} & ... & w_{KD} \\
\end{pmatrix}=\begin{pmatrix}
p_1(x_1) & ... & p_k(x_1) \\
... & ... & ... \\
p_1(x_N) & ... & p_k(x_N) \\
\end{pmatrix}(這裡參看上面的y(x)) = \begin{pmatrix}
p(c_1|x_1) & ... & p(c_k|x_1) \\
... & ... & ... \\
p(c_1|x_N) & ... & p(c_k|x_N) \\
\end{pmatrix}(寫的正規點),\color{red}{這步很重要}\)
\(\tilde X \tilde W - T = \begin{pmatrix} p(c_1|x_1)-t_{11} & ... & p(c_k|x_1)-t_{1K} \\ ... & ... & ... \\ p(c_1|x_N)-t_{1N} & ... & p(c_k|x_N) -t_{NK}\\ \end{pmatrix}\)
\(令\tilde X \tilde W - T = A,a_{ij}=p(c_j|c_i) - t_{ij}\)
\((\tilde X \tilde W - T)^T(\tilde X \tilde W - T)=\begin{pmatrix}
a_{11} & ... & a_{N1} \\
... & ... & ... \\
a_{1K} & ... & a_{KN}\\
\end{pmatrix}\begin{pmatrix}
a_{11} & ... & a_{1K} \\
... & ... & ... \\
a_{N1} & ... & a_{KN}\\
\end{pmatrix}=\begin{pmatrix}
a_{11}^2 + a_{21}^2+...+a_{N1}^2 & ... & ... \\
... & ... & ... \\
... & ... & a_{1K}^2 + a_{2K}^2+...+a_{NK}^2\\
\end{pmatrix}_{K\times K}\)
\(TR(A^TA)=(a_{11}^2 + a_{21}^2+...+a_{N1}^2) + ... + (a_{1K}^2 + a_{2K}^2+...+a_{NK}^2)=C_1分類誤差之和+...+C_K分類誤差之和\)
證畢
令關於\(\widetilde{W}\)的導數等於零，整理，可得\(\widetilde{W}\)的解
\[\widetilde{W} = (\tilde{X}^T\tilde{X})^{-1}\tilde{X}^TT=\tilde{X}^+T \tag{4.16} \]其中\(\tilde{X}^+\)是3.1.1節中討論過的\(\tilde{X}\)的偽逆。然後我們得到形式為
\[\color{red}{y(x) = \widetilde{W}^T\tilde{x} = T^T\left(\tilde{X}^+\right)^T\tilde{x} }\tag{4.17} \]的判別函式。
多目標變數的最小二乘解的一個有趣性質是，如果訓練集合中的每個目標向量對於某些常數\(a, b\)都滿足線性約束
$ a^Tt_n + b = 0 \tag{4.18} $
那麼，模型對於任意的\(x\)值預測也滿足這個約束，即
$ a^Ty(x) + b = 0 \tag{4.19} $
證明 4.18 習題4.2
python 最小二乘法分類
點選檢視程式碼
# 二分類
x_train, y_train = create_toy_data()
x1_test, x2_test = np.meshgrid(np.linspace(-5, 5, 100), np.linspace(-5, 5, 100))
x_test = np.array([x1_test, x2_test]).reshape(2, -1).T
feature = PolynomialFeature(1)
X_train = feature.transform(x_train)
X_test = feature.transform(x_test)
model = LeastSquaresClassifier()
model.fit(X_train, y_train)
y = model.classify(X_test)
plt.scatter(x_train[:, 0], x_train[:, 1], c=y_train)
plt.contourf(x1_test, x2_test, y.reshape(100, 100), alpha=0.2, levels=np.linspace(0, 1, 3))
plt.xlim(-5, 5)
plt.ylim(-5, 5)
plt.gca().set_aspect('equal', adjustable='box')
plt.show()
點選檢視程式碼
import numpy as np
from prml.linear.classifier import Classifier
from prml.preprocess.label_transformer import LabelTransformer
class LeastSquaresClassifier(Classifier):
"""
Least squares classifier model
X : (N, D)
W : (D, K)
y = argmax_k X @ W
"""
def __init__(self, W:np.ndarray=None):
self.W = W
def fit(self, X:np.ndarray, t:np.ndarray):
"""
least squares fitting for classification
Parameters
----------
X : (N, D) np.ndarray
training independent variable
t : (N,) or (N, K) np.ndarray
training dependent variable
in class index (N,) or one-of-k coding (N,K)
"""
if t.ndim == 1:
t = LabelTransformer().encode(t)
self.W = np.linalg.pinv(X) @ t
def classify(self, X:np.ndarray):
"""
classify input data
Parameters
----------
X : (N, D) np.ndarray
independent variable to be classified
Returns
-------
(N,) np.ndarray
class index for each input
"""
return np.argmax(X @ self.W, axis=-1)
點選檢視程式碼
import numpy as np
class LabelTransformer(object):
"""
Label encoder decoder
Attributes
----------
n_classes : int
number of classes, K
"""
def __init__(self, n_classes:int=None):
self.n_classes = n_classes
@property
def n_classes(self):
return self.__n_classes
@n_classes.setter
def n_classes(self, K):
self.__n_classes = K
self.__encoder = None if K is None else np.eye(K)
@property
def encoder(self):
return self.__encoder
def encode(self, class_indices:np.ndarray):
"""
encode class index into one-of-k code
Parameters
----------
class_indices : (N,) np.ndarray
non-negative class index
elements must be integer in [0, n_classes)
Returns
-------
(N, K) np.ndarray
one-of-k encoding of input
"""
if self.n_classes is None:
self.n_classes = np.max(class_indices) + 1
return self.encoder[class_indices]
def decode(self, onehot:np.ndarray):
"""
decode one-of-k code into class index
Parameters
----------
onehot : (N, K) np.ndarray
one-of-k code
Returns
-------
(N,) np.ndarray
class index
"""
return np.argmax(onehot, axis=1)
點選檢視程式碼
x_train, y_train = create_toy_data(add_outliers=True)
x1_test, x2_test = np.meshgrid(np.linspace(-5, 15, 100), np.linspace(-5, 15, 100))
x_test = np.array([x1_test, x2_test]).reshape(2, -1).T
feature = PolynomialFeature(1)
X_train = feature.transform(x_train)
X_test = feature.transform(x_test)
least_squares = LeastSquaresClassifier()
least_squares.fit(X_train, y_train)
y_ls = least_squares.classify(X_test)
logistic_regression = LogisticRegression()
logistic_regression.fit(X_train, y_train)
y_lr = logistic_regression.classify(X_test)
plt.subplot(1, 2, 1)
plt.scatter(x_train[:, 0], x_train[:, 1], c=y_train)
plt.contourf(x1_test, x2_test, y_ls.reshape(100, 100), alpha=0.2, levels=np.linspace(0, 1, 3))
plt.xlim(-5, 15)
plt.ylim(-5, 15)
plt.gca().set_aspect('equal', adjustable='box')
plt.title("Least Squares")
plt.subplot(1, 2, 2)
plt.scatter(x_train[:, 0], x_train[:, 1], c=y_train)
plt.contourf(x1_test, x2_test, y_lr.reshape(100, 100), alpha=0.2, levels=np.linspace(0, 1, 3))
plt.xlim(-5, 15)
plt.ylim(-5, 15)
plt.gca().set_aspect('equal', adjustable='box')
plt.title("Logistic Regression")
plt.show()
python 最小二乘法三分類
點選檢視程式碼
# 三分類
x_train, y_train = create_toy_data(add_class=True)
x1_test, x2_test = np.meshgrid(np.linspace(-5, 10, 100), np.linspace(-5, 10, 100))
x_test = np.array([x1_test, x2_test]).reshape(2, -1).T
feature = PolynomialFeature(1)
X_train = feature.transform(x_train)
X_test = feature.transform(x_test)
least_squares = LeastSquaresClassifier()
least_squares.fit(X_train, y_train)
y_ls = least_squares.classify(X_test)
logistic_regression = SoftmaxRegression()
logistic_regression.fit(X_train, y_train, max_iter=1000, learning_rate=0.01)
y_lr = logistic_regression.classify(X_test)
plt.subplot(1, 2, 1)
plt.scatter(x_train[:, 0], x_train[:, 1], c=y_train)
plt.contourf(x1_test, x2_test, y_ls.reshape(100, 100), alpha=0.2, levels=np.array([0., 0.5, 1.5, 2.]))
plt.xlim(-5, 10)
plt.ylim(-5, 10)
plt.gca().set_aspect('equal', adjustable='box')
plt.title("Least squares")
plt.subplot(1, 2, 2)
plt.scatter(x_train[:, 0], x_train[:, 1], c=y_train)
plt.contourf(x1_test, x2_test, y_lr.reshape(100, 100), alpha=0.2, levels=np.array([0., 0.5, 1.5, 2.]))
plt.xlim(-5, 10)
plt.ylim(-5, 10)
plt.gca().set_aspect('equal', adjustable='box')
plt.title("Softmax Regression")
plt.show()
5-Fisher線性判別函式
我們可以從降維的角度來觀察線性分類模型。首先考慮二分類的情形，並假設有\(D\)維輸入向量\(x\)並使用
$ y=w^Tx \tag{4.20} $
把它投影到一維空間。如果我們在\(y\)上放一個閾值，並把\(y \geq -w_0\)分到類別\(C_1\)中，其餘的分到類別\(C_2\)，那麼我們就得到了前一小節討論的標準的線性分類器。通常來說，投影在一維空間會造成大程度的資訊損失，所以能夠在原來的\(D\)維空間中完全分開的樣本在一維空間中可能會相互重疊。不過通過調節分量的權向量\(w\)，我們可以選擇使類別最分散的一個投影。首先，考慮一個包含\(N_1\)個\(C_1\)點和\(N_2\)個\(C_2\)點的二分類問題，那麼兩個類別的均值向量為：
\[\textbf{m}_1 = \frac{1}{N_1}\sum\limits_{n \in C_1}x_n , \textbf{m}_2 = \frac{1}{N_2}\sum\limits_{n \in C_2}x_n \tag{4.21} \]當投影到\(w\)上後，最簡單的度量類別之間分開程度的方法就是投影之後的的類別均值的距離。這要求我們去選擇那個使得
$ m_2-m_1 = w^T(\textbf{m}_2 - \textbf{m}_1) \tag{4.22} $
\(\textbf{m}_2 , \textbf{m}_1投影在y=w^Tx上的點的距離相減\)
最大的\(w\)的值。 其中
$ m_k = w^T\textbf{m}_k \tag{4.23} $
類別\(C_k\)的資料投影后的均值。但是，可以通過增大\(w\)來使這個表示式無窮大。為了解決這個問題，我們現在\(w\)為單位長的，即\(\sum_i w_i^2 = 1\)。使用拉格朗日乘數法來求解限制條件下的最大化問題，得到\(w \propto (m_2 - m_1)\)。但是這個方法還有圖4.6展示的問題。
證明 \(w \propto (m_2 - m_1)\) ,習題4.4
\(有m_2 -m_1= w^T(\textbf{m}_2 - \textbf{m}_1)\)
\(有w^Tw=1,求w \propto (m_2 - m_1)\)
\(L(w,\lambda)=w^T(m2-m1)+\lambda(w^Tw-1)\)
\(\frac{\partial L }{\partial w}=m_2-m_1 +2\lambda w=0\)
\(w=-\frac{1}{2\lambda}(m_2-m_1)\)
\(w \propto (m_2 - m_1)\)
圖4.6,它展示了在原來二維空間\((x1, x2)\)中可以完全分開的兩個類別，當投影到連線它們的均值的直線上時，有一定程度的重疊。這是由於類別分佈的協方差強非對角引起的(概率分佈的協方差矩陣與對角化矩陣差距較大（暫且沒有推導）)。Fisher提出了最大化一種在給出大的投影均值的距離的同時，給出較小的每個類別內部的方差的函式的主意，來最小化類別重疊。（類內小，類間大）
投影公式（4.20）把通過\(x\)標記的資料點轉換為通過一維空間\(y\)標記的點。轉換後的類別\(C_k\)的內部方差由
$ s_k^2 = \sum\limits_{n \in C_k}(y_n - m_k)^2 \tag{4.24} $
其中\(y_n = w^Tx_n\)。我們可以簡單的定義整個資料集總的類別內部的方差為\(s_1^2 + s_2^2\)。Fisher準則是由類別間的方差與類別內部的方差比例定義的，即
\[J(w) = \frac{(m_2 - m_1)^2}{s_1^2 + s_2^2} \tag{4.25} \]使用式（4.20）、式（4.23）和式（4.24）對它重寫得到
\[J(w) = \frac{w^TS_Bw}{w^TS_Ww} \tag{4.26} \]來顯示的表達對\(w\)的依賴。其中\(S_B\)是類別間的協方差矩陣。\(\color{red}{S_B 類間方差，B代表Between}\)，由
\[S_B = (m_2 - m_1)(m_2 - m_1)^T \tag{4.27} \]給出，\(S_W\)是總得類別內部方差矩陣,\(\color{red}{S_W,類內方差，W代表 within}\)，由
$ S_W = \sum\limits_{n \in C_1}(x_n - m_1)(x_n - m_1)^T + \sum\limits_{n \in C_2}(x_n - m_2)(x_n - m_2)^T \tag{4.28} $
證明 4.26式 習題 4.5
4.20
\(y=w^Tx\)
4.23
\(m_k = w^T\textbf{m}_k\)
4.24
\(s_k^2 = \sum\limits_{n \in C_k}(y_n - m_k)^2\)
\(m_1=w^T\textbf{m}_1\)
\(m_2=w^T\textbf{m}_2\)
\(s_1^2 = \sum\limits_{n \in C_1}(w^Tx_n - w^Tm_1)^2\)
\(s_2^2 = \sum\limits_{n \in C_2}(w^Tx_n - w^Tm_2)^2\)
\(J(w)=\frac{(w^T(m_2-m_1))((m_2-m_1)^Tw)}{\sum\limits_{n \in C_1}(w^T(x_n-m_1))((x_n-m_1)^Tw) \sum\limits_{n \in C_2}(w^T(x_n-m_2))((x_n-m_2)^Tw)}\)
\(S_B=(m_2-m_1)(m_2-m_1)^T\)
\(S_W=\sum\limits_{n \in C_1}((x_n-m_1))((x_n-m_1)^T) \sum\limits_{n \in C_2}((x_n-m_2))((x_n-m_2)^T)\)
對式（4.26）對於\(w\)求導，得到當
$ (w^TS_Bw)S_Ww = (w^TS_Ww)S_Bw \tag{4.29} $
使得\(J(w)\)最大。從式（4.27）我們得到\(S_Bw\)總是在\((m_2 - m_1)\)的方向上。更重要的是我們只在乎\(w\)的方向，而不在乎它的大小，因此我們可以去掉標量因子\((w^TS_Bw) (w^TS_Ww)\) 。在式（4.29）兩邊都乘以\(S_W^{-1}\)得到
$ w \propto S_W^{-1}(m_2 - m_1) \tag{4.30} $
注意，如果類別內部的協方差是各向同性的，那麼\(S_W\)正比於單位矩陣，然後得到和之前討論的一樣，\(w\)正比於類均值的差。
儘管嚴格來說式（4.30）並不是一個判別式，而是對於資料向一維投影的方向的一個具體選擇（得到方向w），但是仍習慣把它稱為Fisher線性判別式。然而，投影的資料可以通過選擇一個閾值\(y_0\)，使得當\(y(x) \geq y_0\)時，我們把它分到\(C_1\)中，否則把它分到\(C_2\)，來構造後面的判別式。
例如，我們可以使用高斯概率分佈對類條件密度$ p(y|C_k) $建模，然後使用1.2.4節(高斯分佈)的技術通過最大似然找到高斯分佈的引數值。
當找到投影類別的高斯近似之後，1.5.1節的方法給出了最優的閾值的形式化解（最小化錯誤分類率）。注意，$ y = w^Tx $是一組隨機變數的和，因此根據中心極限定理，我們可以作出高斯分佈的假設。
這裡提下,在多元統計分析中，獲得了方向\(w\)就可以了，偏置\(w_0\)無所謂，見https://www.cnblogs.com/boyknight/p/16047144.html
4.26求導
借用公式
\(\frac{\partial}{\partial x}\frac{(Ax)^TAx}{(Bx)^TBx}=2\frac{A^TAx}{(Bx)^TBx}-2\frac{(x^TA^TAx)B^TBx}{(x^TB^TBx)^2}--這個公式不知道哪裡得到的，哎，還是要一些數學的基礎功底的\)
\(\frac{\partial}{\partial w}J(w)=2\frac{S_Bw}{w^TS_ww}-2\frac{(w^TS_Bw)S_ww}{(w^TS_ww)^2}\)
\((w^TS_Bw)S_ww =(w^TS_ww)S_Bw\)
Fisher判別 python實現
點選檢視程式碼
# Fisher
x_train, y_train = create_toy_data()
x1_test, x2_test = np.meshgrid(np.linspace(-5, 5, 100), np.linspace(-5, 5, 100))
x_test = np.array([x1_test, x2_test]).reshape(2, -1).T
model = FishersLinearDiscriminant()
model.fit(x_train, y_train)
y = model.classify(x_test)
plt.scatter(x_train[:, 0], x_train[:, 1], c=y_train)
plt.contourf(x1_test, x2_test, y.reshape(100, 100), alpha=0.2, levels=np.linspace(0, 1, 3))
plt.xlim(-5, 5)
plt.ylim(-5, 5)
plt.gca().set_aspect('equal', adjustable='box')
plt.show()
點選檢視程式碼
import numpy as np
from prml.linear.classifier import Classifier
from prml.rv.gaussian import Gaussian
class FishersLinearDiscriminant(Classifier):
"""
Fisher's Linear discriminant model
"""
def __init__(self, w:np.ndarray=None, threshold:float=None):
self.w = w
self.threshold = threshold
def fit(self, X:np.ndarray, t:np.ndarray):
"""
estimate parameter given training dataset
Parameters
----------
X : (N, D) np.ndarray
training dataset independent variable
t : (N,) np.ndarray
training dataset dependent variable
binary 0 or 1
"""
X0 = X[t == 0]
X1 = X[t == 1]
m0 = np.mean(X0, axis=0)
m1 = np.mean(X1, axis=0)
cov_inclass = np.cov(X0, rowvar=False) + np.cov(X1, rowvar=False)
self.w = np.linalg.solve(cov_inclass, m1 - m0)
self.w /= np.linalg.norm(self.w).clip(min=1e-10)
g0 = Gaussian()
g0.fit((X0 @ self.w))
g1 = Gaussian()
g1.fit((X1 @ self.w))
root = np.roots([
g1.var - g0.var,
2 * (g0.var * g1.mu - g1.var * g0.mu),
g1.var * g0.mu ** 2 - g0.var * g1.mu ** 2
- g1.var * g0.var * np.log(g1.var / g0.var)
])
if g0.mu < root[0] < g1.mu or g1.mu < root[0] < g0.mu:
self.threshold = root[0]
else:
self.threshold = root[1]
def transform(self, X:np.ndarray):
"""
project data
Parameters
----------
X : (N, D) np.ndarray
independent variable
Returns
-------
y : (N,) np.ndarray
projected data
"""
return X @ self.w
def classify(self, X:np.ndarray):
"""
classify input data
Parameters
----------
X : (N, D) np.ndarray
independent variable to be classified
Returns
-------
(N,) np.ndarray
binary class for each input
"""
return (X @ self.w > self.threshold).astype(np.int)
6-最小平方和Fisher判別的關係
確定線性判別式的最小二乘方法是基於使模型預測儘可能的接近目標值的目的的。相反，Fisher準則的目標是最大化輸出空間中類別的區分度。這兩種方法之間的關係是很有趣的。特別的，我們會證明，在二分類問題中，Fisher準則可以看成最小二乘的一個特例。
目前為止，我們一直採用“1-of-K”編碼來表示目標值。然而，如果我們採用一種稍微不同的編碼方式，那麼權重的最小二乘解會等價於Fisher判別式的解（Duda and Hart, 1973）。特別的，我們讓屬於\(C_1\)的目標值等於\(N/N_1\)，其中\(N_1\)是類別\(C_1\)的模式的數量，\(N\)是總的模式數量。這個目標值近似於類別\(C_1\)的先驗概率的倒數，同時令\(C_2\)目標值等於\(−N/N_2\)，其中\(N_2\)是類別\(C_2\)的模式的數量。
平方和誤差函式可以寫成
分別關於\(w_0, w\)求\(E\)的導數，並使其等於0，得到
\[\begin{eqnarray} \sum\limits_{n=1}^N\left(w^Tx_n + w_0 - t_n \right) = 0 \tag{4.32} \\ \sum\limits_{n=1}^N\left(w^Tx_n + w_0 - t_n \right)x_n = 0 \tag{4.33} \end{eqnarray} \]根據式（4.32），並按選擇的目標編碼方式來編碼\(t_n\)，就可得到偏置的表示式
$ w_0 = -w^Tm \tag{4.34} $
其中我們使用了
$ \sum\limits_{n=1}^Nt_n = N_1\frac{N}{N_1} - N_2\frac{N}{N_2} = 0 \tag{4.35} $
且\(m\)是由
$ m = \frac{1}{N}\sum\limits_{n=1}^Nx_n = \frac{1}{N}(N_1m_1 + N_2m_2) \tag{4.36} $
給出的全部資料的均值。
說白了就是兩個演算法的編碼方式不一樣