[ML筆記]Logistic Regression

scikit-learn

使用 scikit-learn 訓練一個感知器

from sklearn import datasets
import numpy as np

iris = datasets.load_iris() # 就是那堆花
X = iris.data[:, [2, 3]] # 選取花瓣長度(2)與花瓣寬度(3)當作特徵
y = iris.target # 花的類別

print('X:', X[1:5, ]) # 檢查print('Label:', np.unique(y))

終端輸出:

> X: [[1.4 0.2]
 [1.3 0.2]
 [1.5 0.2]
 [1.4 0.2]]
> Label: [0 1 2]

在這裡的 X 代表了那堆花的特徵(features),這裡用的是花瓣長度([2])與花瓣寬度([3])當作分類時的特徵,而 y 則是花的類別(本來是 sentosa、versicolor、virginia,這裡用數字代替)。

接著,我們把所有的資料分割成訓練用(train set)的與測試用(test set)的兩坨。使用的是 scikit-learn 裡的 train_test_split() 函式。這裡用的比例是 train set : test set = 3:7,當然要改成別的比例也可以,反正就是分成兩邊。

from sklearn.model_selection import train_test_split
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=1, stratify=y)

接著將資料標準化,用上 scikit-learn 裡的 StandardScaler() 函式。

from sklearn.preprocessing import StandardScaler
sc = StandardScaler()
sc.fit(X_train)
X_train_std = sc.transform(X_train)
X_test_std = sc.transform(X_test)

然後就可以開始訓練感知器惹。這個 scikit-learn 裡的 Perceptron() 函式用法和之前所說的感知器很接近,有一個學習速率(eta0)和迭代次數(n_iter)。

from sklearn.linear_model import Perceptron
ppn = Perceptron(n_iter = 40, eta0 = 0.1)
ppn.fit(X_train_std, y_train)

這樣就分類完了,接下來檢查有幾個東西是「分類錯誤」的。

y_pred = ppn.predict(X_test_std)
print('分類錯誤的樣本: %d' % (y_test != y_pred).sum())
> 分類錯誤的樣本: 9

還有分類準確率(accuracy rate):

from sklearn.metrics import accuracy_score
print('準確率: %.2f' % accuracy_score(y_test, y_pred))
> 準確率: 0.80

這個是這樣算出來的:

\begin{equation} 1 - \textrm{分類錯誤率} = 1 - \frac{9\textrm{(分類錯誤的樣本)}}{150\textrm{(總樣本數)} \times 0.3\textrm{(分割比例)}}=1-0.2=0.8=80\% \end{equation}

視覺化

import matplotlib.pyplot as plt
from matplotlib.colors import ListedColormap

def plot_decision_regions(X, y, classifier, test_idx=None, resolution = 0.02):
	# 設定畫布
    markers = ('s', 'x', 'o', '^', 'v')
    colors = ('red', 'blue', 'lightgreen', 'gray', 'cyan')
    cmap = ListedColormap(colors[:len(np.unique(y))])

    x1_min, x1_max = X[:, 0].min() - 1, X[:, 0].max() + 1
    x2_min, x2_max = X[:, 1].min() - 1, X[:, 1].max() + 1
    xx1, xx2 = np.meshgrid(np.arange(x1_min, x1_max, resolution),
                           np.arange(x2_min, x2_max, resolution))
    Z = classifier.predict(np.array([xx1.ravel(), xx2.ravel()]).T)
    Z = Z.reshape(xx1.shape)
    plt.contourf(xx1, xx2, Z, alpha=0.3, cmap = cmap)
    plt.xlim(xx1.min(), xx1.max())
    plt.ylim(xx2.min(), xx2.max())

    for idx, cl in enumerate(np.unique(y)):
        plt.scatter(x = X[y == cl, 0], 
                    y = X[y == cl, 1],
                    alpha = 0.8, 
                    c = colors[idx],
                    marker = markers[idx], 
                    label = cl, 
                    edgecolor ='black')
	
	# 特別標出要預測的樣本,把他們用圈圈表示
    if test_idx:
        X_test, y_test = X[test_idx, :], y[test_idx]
        plt.scatter(X_test[:, 0],
                    X_test[:, 1],
                    c='',
                    edgecolor ='black',
                    alpha = 1.0,
                    linewidth = 1,
                    marker = 'o',
                    s = 100, 
                    label = 'test set')

X_combined_std = np.vstack((X_train_std, X_test_std))
y_combined = np.hstack((y_train, y_test))

plot_decision_regions(X = X_combined_std, y = y_combined, classifier = ppn, test_idx = range(105, 150))
plt.xlabel('標準化後的花瓣長度')
plt.ylabel('標準化後的花瓣長度')
plt.legend(loc = 'upper left')

plt.tight_layout()
plt.show()

為什會有分類錯誤的漏網之魚呢?那是因為這個感知器最好是只應用於線性可分(linearly separable)的樣本的,對於線性不可分(linearly unseparable)的資料,自然會有沒辦法完美劃分的可能。畢竟線性感知器是最簡單的分類方法,所以實際上很少使用。

Logistic Regression

這其實是一個取名的藝術,Logistic Regression 雖然美其名有的 Regression(回歸),但是他基本上更像是一種分類(Classification)的方法。Logistic Regression 的基本概念是把東西分成 $$1$$ 和 $$0$$ ,對所有的輸入而言,都只回得到兩種可能的輸出情況。

舉例而言,「或許」我們可以透過身高和體重之間的關係(或許啦,我不確定),來判斷這個人是男生還是女生(二元變數,binary variable)。對於只有「身高」和「體重」兩種變數,如果我們把這散佈在二維平面圖上的關係點求出一條線性回歸(linear regression)線,然後用這條回歸線說未知的樣本是男還是女,那就可能出現這個人是「0.7 個男生 + 0.3個女生」這種事情。問題是,我們分類的目標就只有男生跟女生啊。

又或者是我嘗試從今天的天氣判斷出巷口的麵店會不會開,實際上只會有「開」跟「不開」兩種可能性。其他像是「死了」跟「活著」、「生病了」跟「很健康」、「上漲」跟「下跌」,這些分類都是二元對立的。所以我們導入機率的概念,把迴歸係數的解讀為「當自變項增加一個單位,依變項 A 相對依變項 B C 的機率以及依變項 A B 相對依變項 C 的機率會增加幾倍」。

數學推導

我覺得都還給計量老師惹…從條件機率開始吧,今天令條件機率等於一個 $F(\cdot)$ 函數:

\begin{equation} \underbrace{P(Y=1 \vert X ) =\frac{P(Y=1 \cap X)}{P(X)}}_{\textrm{條件機率定義}} \equiv F(w_{0}x_{0}+w_{1}x_{1}+\cdots+w_{n}x_{n}) =F(\textbf{w}^\mathrm{T}\cdot\textbf{x})\\ \textrm{其中,我們用兩個向量表示的話,分別令 }\mathbf{w} = \left(\begin{array}{ccc} w_{0}\\w_{1}\\ \vdots \\ w_{n}\end{array}\right) , \textbf{x} = \left(\begin{array}{ccc}x_{0}\\x_{1}\\ \vdots \\ x_{n}\end{array}\right) \end{equation}

上面條件機率(conditional probability)的定義就是「在給定 $X$ 發生的情況之下,發生 $Y=1$ 的機率是多少」。

而右式的 $F(\cdot)$ 我們定義其為:

\begin{equation} F(\cdot) \equiv F(w_{0}x_{0}+w_{1}x_{1}+\cdots+w_{n}x_{n})=F(\textbf{w}^\mathrm{T}\cdot\textbf{x}) =\frac{1}{1+e^{-(w_{0}x_{0}+w_{1}x_{1}+\cdots+w_{n}x_{n})}} =\frac{1}{1+e^{(-\mathbf{w}^\mathrm{T}\cdot\textbf{x})}} \end{equation}

同理,如果把上面 $F(\cdot)$ 的函數輸入值改以一個淨輸入代替時,我們稱其為 Sigmoid 函數,Sigmoid 函數就是 Logistic Regression 的激勵函數:

\begin{equation} \textrm{當我們令 }z = \textbf{w}^\mathrm{T}\cdot\textbf{x}=w_{0}x_{0}+w_{1}x_{1}+\cdots+w_{n}x_{n} \textrm{時,} \phi(z)=F(z)=\frac{1}{1+e^{-z}}\cdots \textrm{Sigmoid Function} \end{equation}

因為我們的結果只有兩種(剛剛提到的二元對立),所以發生 $Y=0$ 的機率就是一減去 $Y=1$ 發生的機率:

\begin{equation} P(Y=0\vert X)=1-P(Y=1\vert X) \end{equation}

看到這裡,是不是覺得這個東西很伯努利分佈咧,都是:

\begin{equation} \left\{ \begin{array}{l} \textrm{發生 Y=1,或稱 target 時的機率} &= P(Y=1|X)\\ \textrm{發生 Y=0,或稱 non-target 時的機率} &= P(Y=0|X) = 1-P(Y=1|X)\ \end{array} \right . \end{equation}

所以使用最大慨似估計(maximum likelihood estimation, MLE)進行估計:

\begin{equation} \begin{split} \mathcal{L}(\mathbf{w})&=(F(z)^{y_{1}}(1-F(z))^{(1-y_{1})})\times(F(z)^{y_{2}}(1-F(z))^{(1-y_{2})})\times\cdots\times(F(z)^{y_{n}}(1-F(z))^{(1-y_{n})})\\ &=\prod_{i=1}^{n} F(z)^{y_{i}}(1-F(z))^{1-y_{i}}\\ \end{split} \end{equation}

取自然對數 ln 後得到:

\begin{equation} \begin{split} ln\mathcal{L}(\mathbf{w})\equiv \mathcal{l}(\mathbf{w}) &=\ln \Big(\prod_{i=1}^{n} F(z)^{y_{i}}(1-F(z))^{1-y_{i}}\Big)\\ &=\sum_{i=1}^{n}\ln \Big(F(z)^{y_{i}}(1-F(z))^{1-y_{i}}\Big)\\ &=\sum_{i=1}^{n}\Big(\ln(F(z)^{y_{i}})+\ln((1-F(z))^{1-y_{i}}\Big)\\ &=\sum_{i=1}^{n}\Big(y_{i}\ln(F(z))+(1-y_{i})\ln(1-F(z)\Big)\\ &=\sum_{i=1}^{n}\Big(y_{i}\ln(F(z))+\ln(1-F(z)-y_{i}\ln(1-F(z)\Big)\cdots(*)\\ &=\sum_{i=1}^{n}\Big(y_{i}\ln\frac{F(z)}{1-F(z)}+\ln(1-F(z)\Big)\\ \end{split} \end{equation}

接著,取導求極值,要 chain rule:

\begin{equation} \begin{split} \frac{\partial\ln\mathcal{L}(\mathbf{w})}{\partial w_{i}} &=\frac{\partial\mathcal{l}(\mathbf{w})}{\partial w_{i}} =\Big(y_{i}\frac{1}{F(z)}-(1-y_{i})\frac{1}{1-F(z)}\Big)\frac{\partial F(z)}{\partial w_{i}}\\ \textrm{其中,}&\frac{\partial F(z)}{\partial w_{i}}=\frac{\partial z}{\partial w_{i}}\frac{\partial}{\partial z}\Big(\frac{1}{1+e^{-z}}\Big)\\ &=\frac{e^{-z}}{(1+e^{-z})^2}\frac{\partial z}{\partial w_{i}} =\Big(\frac{1}{1+e^{-z}}\Big)\Big(1-\frac{1}{1+e^{-z}}\Big)\frac{\partial z}{\partial w_{i}}\\& =F(z)(1-F(z))\frac{\partial z}{\partial w_{i}}\\ \textrm{代回,}&\Big(y_{i}\frac{1}{F(z)}-(1-y_{i})\frac{1}{1-F(z)}\Big)\frac{\partial F(z)}{\partial w_{i}}\\ &=\Big(y_{i}\frac{1}{F(z)}-(1-y_{i})\frac{1}{1-F(z)}\Big)\Big(F(z)(1-F(z)\Big)\frac{\partial z}{\partial w_{i}}\\ &=\Big(y_{i}(1-F(z))-(1-y_{i}F(z)\Big)x_{i}\\ &=\Big(y_{i}-F(z)\Big)x_{i} \end{split} \end{equation}

我們的目標是讓最大慨似估計最大化的權重值,所以在得到導數後可以加總權重。

\begin{equation} w_{i}:=w_{i}+\eta \underbrace {\sum_{i=1}^{n}(y_{i}-F(z))x_i}_{\equiv \Delta w_{i}} \end{equation}

也可以一次更新所有的權重,得到

\begin{equation} \mathbf{w}:=\mathbf{w}+\Delta\mathbf{w}=\mathbf{w}+\eta\nabla\mathcal{l}(\mathbf{w}) \end{equation}

剛好等價於梯度下降的權重更新公式!

而代價函數則可以用 $-(*)$ 來表示:

\begin{equation} J(\mathbf{w})=-\sum_{i}(y_{i}\ln(F(z))+\ln((1-F(z))-y_{i}\ln((1-F(z))) \end{equation}

畫出來!

嘗試畫一個 Sigmoid 函數。

import matplotlib.pyplot as plt
import numpy as np

def sigmoid(z):
    return 1.0 / (1.0 + np.exp(-z))

z = np.arange(-7, 7, 0.1)
phi_z = sigmoid(z)

plt.plot(z, phi_z)
plt.axvline(0.0, color = 'k')
plt.ylim(-0.1, 1.1)
plt.xlabel('z')
plt.ylabel('$\phi (z)$')

# y 軸標出 y = 0, 0.5, 1 的線
plt.yticks([0.0, 0.5, 1.0])
ax = plt.gca()
ax.yaxis.grid(True)

plt.tight_layout()
plt.show()

訓練

這裡直接使用 sklearn.linear_model 內建的 LogisticRegression 函式。

from sklearn.linear_model import LogisticRegression

lr = LogisticRegression(C = 100.0, random_state = 1)
lr.fit(X_train_std, y_train)

plot_decision_regions(X_combined_std, y_combined,
                      classifier = lr, test_idx = range(105, 150))
plt.xlabel('標準化後的花瓣長度')
plt.ylabel('標準化後的花瓣長度')
plt.legend(loc = 'upper left')
plt.tight_layout()
plt.show()