廣義線性模型與邏輯回歸

02-16

零、概述

本文路線圖：線性模型 -> 廣義線性模型 -> 邏輯回歸 -> 牛頓法以及 IRLS 。前導知識包括：矩陣、向量的計算，多元函數的梯度、赫森矩陣以及二階泰勒展開，概率論。

一、線性模型

假設樣本由 $p$ 維特徵向量 $X$ 和標量響應 $y$ 組成。 $X$ 和 $y$ 服從聯合分布 $Prleft(y, X ight)$ 。把確定 $X=x$ 的前提下 $y$ 的條件分布記作 $Prleft(y|X=x ight)$ 。如果假設對於 $X=x$ ， $y$ 的隨機雜訊是加性的，如下式：

$y = fleft( x ight) + varepsilon left[ 1.1 ight]$

則當 $X=x$ 時 $y$ 本質上由 $fleft(x ight)$ 確定，其所有不確定性來源於 $varepsilon$ 。如果進一步假設 $varepsilon$ 滿足高斯分布：

$varepsilonsim Nleft(0, sigma^2 ight) left[ 1.2 ight]$

那麼這就是加性高斯誤差。在這種設定下分布 $Prleft(y|X=x ight)$ 的條件期望，也是概率密度最大的值就是：

$Eleft(y|X=x ight)=fleft(x ight) left[1.3 ight]$

把該條件期望當作 $X=x$ 時對 $y$ 的預測值。至於 $varepsilon$ 則是問題固有的、無法通過模型優化去除的誤差。所謂線性模型就是用線性函數來建模 $fleft(x ight)$ ：

$Eleft(y|X=x ight)=fleft(x ight) =eta_0+eta_1x_1+eta_2x_2+...+eta_px_p = x^Teta left[1.4 ight]$

最後一個等號後面的 $eta$ 是 $p+1$ 維向量 $left( eta_0, eta_1, ... eta_p ight) ^T$ ， $x$ 是 $p+1$ 維向量 $left(1, x_1, x_2, ... x_p ight)^T$ 。 $x$ 前面加一維常量 1 是為了把截距 $eta_0$ 納入向量表達式中。訓練集（ $training set$ ）是 $N$ 個特徵向量及其標量響應： $left{ x_i, y_i ight}_{i=1}^{N}$ 。如果問題是一個回歸問題，且選用平方誤差作為損失函數：

$Lleft(y,fleft(x ight) ight)=left(y-fleft(x ight) ight)^2 left[1.5 ight]$

令 $hat{eta}$ 在訓練集上最小化平均損失：

$hat{eta} = argmin_etafrac{1}{N}sum_{i=1}^{N}{left(y_i-fleft(x_i ight) ight)^2}=argmin_etafrac{1}{N}sum_{i=1}^{N}{left(y_i-x_i^Teta ight)^2} = argmin_etafrac{1}{N} left( y-Xeta ight)^Tleft( y-Xeta ight) left[1.6 ight]$

最後一個等號後面的 $y$ 是所有 $y_i$ 組成的 $N$ 維向量。 $X$ 是 $N imes p$ 矩陣，其第 $i$ 行是 $x_i^T$ 。損失函數作為 $eta$ 的函數是凸的，存在全局唯一最小點，也就是唯一的駐點。通過求導並置導向量為 0 再解方程可得到：

$hat{eta}=left(X^TX ight)^{-1}X^Ty left[1.7 ight]$

這就是要尋找的最優參數——線性回歸的最小二乘解。當 $eta=hat{eta}$ 時對 $N$ 個訓練樣本的預測值是：

$hat{y}=Xhat{eta}=Xleft(X^TX ight)^{-1}X^Ty =Sy left[1.8 ight]$

$hat{y}$ 是 $X$ 的列的線性組合，係數向量是 $hat{eta}$ 。也就是說 $hat{y}$ 屬於 $X$ 的列空間。並且由於 $hat{eta}$ 最優化了平均平方誤差，所以 $hat{y}$ 是 $X$ 的列空間中與向量 $y$ 歐氏距離最小的向量，即 $hat{y}$ 是 $y$ 在 $X$ 列空間上的投影。同時 $hat{y}$ 也是 $y$ 的一個線性變換，變換矩陣是 $S$ 。這也是該模型被稱為線性模型的一個原因。 $S$ 的對角線元素的和，即 $S$ 的跡為：

$Trleft(S ight)=Trleft(Xleft(X^TX ight)^{-1}X^T ight)=Trleft(left(X^TX ight)^{-1}X^TX ight)=Trleft(I ight)=p left[1.9 ight]$

上式第二個等號利用了 $AB$ 的跡與 $BA$ 的跡相等。把矩陣乘積的跡的表達式按元素展開很容易證明這一點。 $S$ 的跡等於自由參數的個數。所以把 $S$ 的跡定義為線性模型的自由度是合適的。如果為模型加上正則化，那麼 $eta$ 就受到了一種約束，它在參數空間中的變化就沒有那麼「自由」了。正則化的線性模型的自由度就不應該是單純的參數個數。其變換矩陣的跡小於 $p$ ，用跡度量正則化線性模型的自由度仍是恰當的（參見本專欄文章理解機器學習中的 L2 正則化）。

二、廣義線性模型和邏輯回歸

把 $y$ 在 $X=x$ 時的條件期望簡記作：

$muleft(x ight)=Eleft(y|X=x ight) left[ 2.1 ight]$

如果引入一個函數 $g$ ，令：

$gleft(muleft(x ight) ight)=eta_0+eta_1x_1+eta_2x_2+...+eta_px_p = x^Teta left[2.2 ight]$

$g$ 稱作鏈接函數 $Link Function$ 。如果 $g$ 是恆等函數 $gleft(t ight)=t$ ，則 [2.2] 就退化為 [1.4] 。所以線性模型是 [2.2] 在鏈接函數為恆等函數時的特例。[2.2] 就是廣義線性模型。鏈接函數 $g$ 採用不同的形式就產生不同的模型。

考慮二分類問題。樣本輸出 $y$ 是類別的標記：0 表示負例；1 表示正例。把 $X=x$ 時樣本屬於正例的概率記為 $Prleft(y=1|X=x ight)$ ，則有：

$muleft(x ight)=Eleft(Y|X=x ight)=1 imes Prleft(y=1|X=x ight) + 0 imes Prleft(y=0|X=x ight) = Prleft(y=1|X=x ight)=pleft(x ight) left[ 2.3 ight]$

最後的 $pleft(x ight)$ 是 $Prleft(y=1|X=x ight)$ 的簡寫。令鏈接函數 $g$ 是 $logit$ 函數：

$logitleft(t ight) = logfrac{t}{1-t} left[2.4 ight]$

$logit$ 函數圖形為：

圖 2.1 logit 函數圖形

把 $logit$ 函數和 [2.3] 代入 [2.2] 中有：

$logitleft(mu left(x ight) ight)=logfrac{muleft(x ight)}{1-muleft(x ight)}=logfrac{pleft(x ight)}{1-pleft(x ight)}=x^Teta left[2.5 ight]$

也就是將正例與負例概率之比的對數（稱作 $log odds$ ）建模為線性。從這裡也可看出，不論最終將判斷正例的概率閾值定在多少，模型判斷正例與負例的分界線是 $logitleft(frac{pleft(x ight)}{1-pleft(x ight)} ight)=x^Teta=C$ ，為一條直線。現在將 [2.5] 變形得到：

$pleft(x;eta ight)=frac{1}{1+e^{-x^Teta}} left[ 2.6 ight]$

這裡的記法 $pleft(x;eta ight)$ 強調一下 $pleft(x ight)$ 受 $eta$ 影響。[2.6] 就是邏輯回歸模型。可見邏輯回歸屬於廣義線性模型。分類問題與回歸問題表面上看是不同的問題，但在深層次上它們是統一的——都是用損失函數衡量擬合的損失，調整參數盡量降低這個損失。

對於回歸問題，可以用平方誤差 [1.5] 衡量擬合的損失。當模型輸出不等於真實值時應該加以懲罰。平方誤差在模型輸出不等於真實值時大於 0 ，且相差越大則平方誤差越大。但是平方誤差隨著模型預測與真實值差距的增大呈二次增長，對於大差距懲罰得太過。如果存在離群樣本，比如測量錯誤或標註錯誤，平方誤差易受干擾，並不健壯。存在其他懲罰函數以緩和這個問題。分類問題的損失函數又應該有另外的形式。

對於分類問題，現在令 $y$ 取 1 或 -1 標記正負例。模型 $fleft(x ight)$ 是連續函數。如果當 $fleft(x ight)>0$ 時將樣本判為正例；當 $fleft(x ight)leq0$ 時判為負例。那麼當 $y$ 和 $fleft(x ight)$ 同號，即 $ycdot fleft(x ight)>0$ 時模型分類正確；當 $ycdot fleft(x ight)leq 0$ 時模型分類錯誤。注意當 $y$ 只能取 1 或 -1 時，平方誤差可以視作 $ycdot fleft(x ight)$ 的函數：

$Lleft(y,fleft(x ight) ight)=left(y-fleft(x ight) ight)^2=y^2-2ycdot fleft(x ight)+f^2left(x ight) =1-2ycdot fleft(x ight)+left(ycdot fleft(x ight) ight)^2=left(1-ycdot fleft(x ight) ight)^2 left[2.7 ight]$

其圖形為（藍色曲線）：

圖 2.2 三種損失函數的圖形。黑色虛線是錯分類懲罰：分錯懲罰為 1.0；分對懲罰為 0.0 。

觀察藍色曲線。當 $ycdot fleft(x ight)<=0$ 時分類錯誤。如果把 $fleft(x ight)$ 的絕對值作為分類置信程度的話，負的 $y cdot fleft(x ight)$ 越小錯誤越嚴重。藍色線條在這種情況下升高，懲罰這種現象。當 $y cdot fleft(x ight)>0$ 時分類是正確的。且 $y cdot fleft(x ight)$ 越大則分得越對。這時懲罰應該減小。但是平方誤差的行為是先減小，到大於 1.0 後反而增大。平方誤差是評價回歸質量的。回歸問題小也不行大也不行。但是分類問題性質不同，平方誤差不再恰當。

對於二分類問題，一個可採用的損失函數是 $Binomial Deviance$ 。如果 $y$ 取 1 和 0 標記正負例。令 $X=x$ 時樣本為正例的概率是 $pleft(x;eta ight)$ ，該損失函數為：

$Binomial Deviance = -y imes log pleft(x;eta ight)-(1-y) imes log left(1-pleft(x;eta ight) ight) left[2.8 ight]$

它是「對數似然」的相反數。將 [2.8] 取相反就是對數似然。可以到看當 $y=1$ 時對數似然為 $log pleft(x;eta ight)$ ；當 $y=0$ 時對數似然為 $log left(1-pleft(x;eta ight) ight)$ 。 $Binomial Deviance$ 的相反數同時也是交叉熵 $Cross Entropy$ 。它衡量兩個分布的相似程度。其中一個分布是模型的預測分布。另一個分布是觀察到的該樣本類別的分布。對數似然和交叉熵越大越好。最小化 $Binomial Deviance$ 就是最大化對數似然和交叉熵。

如果用 $y=1 or -1$ 來標記正負例，則 $frac{y+1}{2}$ 在正負例時取值分別是 1 和 0 。將 $frac{y+1}{2}$ 和式 [2.6] 帶入 [2.8] ，令 $fleft(x ight)=x^Teta$ ，有：

$Binomial Deviance=-frac{y+1}{2} logfrac{1}{1+e^{-fleft(x ight)}}-frac{1-y}{2}logfrac{1}{1+e^{fleft(x ight)}}=logleft(1+e^{-yfleft(x ight)} ight) left[ 2.9 ight]$

第二步推導有些跳。記住 $y$ 的值只有兩種可能 1 或 -1 。推導不下去時把式中的 $y$ 替換成 1 和 -1 分別繼續推導。最後兩條路合成一個最終的結果就是 [2.9] 。 $Binomial Deviance$ 的圖形是圖 2.2 中的橙色曲線：當 $ycdot fleft(x ight)<=0$ 時增大，但沒有誤差平方增大得那麼快。所以它對於離群樣本更健壯。當 $ycdot fleft(x ight)>0$ 時減小，且隨著 $ycdot fleft(x ight)$ 增大越來越小。這說明對於分類問題來說 $Binomial Deviance$ 的懲罰行為是恰當的。

以下繼續用 $y=1 or 0$ 標記正負例。採用 $Binomial Deviance$ 作為損失函數。邏輯歸回就是尋找 $eta$ 在訓練集上最小化平均損失：

$hat{eta} = argmin_etafrac{1}{N}sum_{i=1}^{N}{left{-y logfrac{1}{1+e^{-x_i^Teta}} -left(1-y ight) logfrac{1}{1+e^{x_i^Teta}} ight}} left[2.10 ight]$

當找到了 $hat{eta}$ ，一個新的樣本 $x$ 屬於正例的概率就是：

$hat{p}left(x;hat{eta} ight)=frac{1}{1+e^{-x^That{eta}}} left[2.11 ight]$

可以根據 $hat{p}$ 是否大於 0.5 （或者其他閾值）判斷 $x$ 是否屬於正例。從這個式子也可以看出邏輯回歸模型相當於只有一個神經元的神經網路，該神經元激活函數是 $Logistic Function$ 。 $Binomial Deviance$ 不是凸函數，其最小值沒有解析解，需要通過最優化演算法來尋找 $hat{eta}$ ，比如梯度下降法（參見本專欄文章神經網路之梯度下降與反向傳播（上））。本文介紹基於目標函數二階性質的牛頓法來尋找邏輯回歸的解。

三、牛頓法與迭代重加權最小二乘（IRLS）

如果一個函數 $Lleft(eta ight)$ 是二階可導的，在任一個數據點 $eta^{old}$ 附近可以將它二階泰勒展開：

$Lleft(eta-eta^{old} ight) approx left(G_L|_{eta^{old}} ight)^Tleft(eta-eta^{old} ight)+frac{1}{2}left(eta-eta^{old} ight)^TH_L|_{eta^{old}}left(eta-eta^{old} ight) left[3.1 ight]$

$Lleft(eta ight)$ 是 $R^p ightarrow R^1$ 函數。 $eta$ 和 $eta^{old}$ 都是 $p$ 維向量。 $G_L|_{eta^{old}}$ 是 $Lleft(eta ight)$ 在 $eta^{old}$ 處的梯度，是一個 $p$ 維向量。 $H_L|_{eta^{old}}$ 是 $Lleft(eta ight)$ 在 $eta^{old}$ 處的二階導矩陣——赫森矩陣，是一個 $p imes p$ 的矩陣。牛頓法的原理就是在一個點將函數二階泰勒展開，展開後得到一個凸二次函數，是原函數在該點附近的二階近似。將該二階近似的全局最小點作為迭代的下一個點。如圖所示：

圖 3.1 牛頓法一步迭代示意圖

二次函數在全局最小點梯度為 0 。將 [3.1] 約等號右側求導後置 0 並解方程求得全局最小點是：

$eta^{new}=eta^{old}-left(H_L|_{eta^{old}} ight)^{-1}G_L|_{eta^{old}} left[3.2 ight]$

[3.2] 就是牛頓法的迭代公式。如果 $Lleft(eta ight)$ 本身就是一個二次函數（例如平方誤差 [1.6] ），則它在任一點的二階泰勒展開就是它本身。這種情況下牛頓法只用一步迭代就可以到達全局最小點。最小二乘解 [1.7] 就是牛頓法一步達到的 $eta^{new}$ 。

現在應用牛頓法來尋找邏輯回歸模型的最優參數。為了方便去掉了 $frac{1}{N}$ 係數，這不影響尋找最小值：

$Binomial Deviance=Lleft(eta ight) \= -sum_{i=1}^{N}left{{y_i log pleft(x_i;eta ight)} + left( 1- y_i ight) log left( 1-pleft( x_i;eta ight) ight) ight} \=-sum_{i=1}^{N} left{{y_i log frac{pleft(x;eta ight)}{1-pleft(x;eta ight)} +log left(1-pleft(x;eta ight) ight)} ight} \=-sum_{i=1}^{N}{left{y_ix_i^Teta-log left(1+e^{x_i^Teta} ight) ight}} left[3.3 ight]$

將 [3.3] 式對 $eta$ 求導（並轉置）得到 $Lleft(eta ight)$ 在 $eta$ 的梯度：

$G_L|_eta=left(frac{delta L}{deltaeta} ight)^T=-sum_{i=1}^{N}{x_ileft(y_i-pleft(x_i;eta ight) ight)} left[3.4 ight]$

將梯度求導得到 $Lleft(eta ight)$ 在 $eta$ 的赫森矩陣：

$H_L|_eta=frac{delta^2 L}{delta eta delta eta^T}=sum_{i=1}^{N}{x_ix_i^T}pleft(x_i;eta ight)left(1-pleft(x_i;eta ight) ight) left[3.5 ight]$

[3.5] 有些跳步，因為中間過程的式子不好編輯。推導時注意 $pleft(x;eta ight)$ 對 $x^Teta$ 的導數是 $pleft(x;eta ight)left(1-pleft(x;eta ight) ight)$ 。[3.4] 和 [3.5] 可以整理成矩陣形式：

$egin{equation} left{ egin{aligned} G_L|_eta=-X^Tleft(y-p ight) \H_L|_eta=X^TWX end{aligned} ight.egin{aligned} quad left[3.6 ight]end{aligned} end{equation}$

$X$ 為 $N imes p$ 矩陣，第 $i$ 行為 $x_i^T$ 。 $y$ 和 $p$ 分別是由 $y_i$ 和 $pleft(x_i;eta ight)$ 組成 $N$ 維向量。 $W$ 是 $N imes N$ 對角矩陣，對角線元素是 $pleft(x_i;eta ight)left(1-pleft(x_i;eta ight) ight)$ 。將 [3.6] 代入牛頓法迭代公式 [3.2] 得到：

$eta^{new}=eta^{old}-left(H_L|_{eta^{old}} ight)^{-1}G_L|_{eta^{old}}=eta^{old}+left(X^TWX ight)^{-1}X^Tleft(y-p ight) left[3.7 ight]$

[3.7] 就是用牛頓法訓練邏輯回歸模型的迭代公式。訓練過程直到 $eta$ 收斂時停止。繼續將 [3.7] 變形：

$eta^{new}=eta^{old}-left(H_L|_{eta^{old}} ight)^{-1}G_L|_{eta^{old}}\=eta^{old}+left(X^TWX ight)^{-1}X^Tleft(y-p ight)\=left(X^TWX ight)^{-1}X^TWXeta^{old}+left(X^TWX ight)^{-1}X^Tleft(y-p ight)\=left(X^TWX ight)^{-1}X^TWleft(Xeta^{old}+W^{-1}left(y-p ight) ight)\=left(X^TWX ight)^{-1}X^TWz left[3.8 ight]$

[3.8] 中的 $z=Xeta^{old}+W^{-1}left(y-p ight)$ 。可以看出 $eta^{new}$ 是以 $X$ 為特徵以 $z$ 為響應的加權最小二乘解。 $z$ 被稱為工作響應。 $z$ 每一次迭代都有變化，是迭代過程中的臨時響應。權值矩陣 $W$ 也隨著每次迭代而變化。故用 [3.8] 求邏輯回歸最優參數又被稱為迭代重加權最小二乘，Iteratively Reweighted Least Squares，IRLS 。

四、實現及舉例

本節採用 kaggle 上關於骨盆疾病的數據。共有六個數值型特徵：pelvic_incidence，pelvic_tilt numeric，lumbar_lordosis_angle，sacral_slope，pelvic_radius，degree_spondylolisthesis 。都是關於骨盆和脊椎的測量值。類別分兩種 Normal 和 Abnormal 。數據下載：kaggle 里的 column_2C_weka.csv 。下列圖表是對數據所做的簡單探查。

圖 4.1 六個特徵分組成對散點圖

圖 4.2 異常與正常樣本的六個特徵的箱狀圖

圖 4.2 異常與正常樣本的數量對比

圖 4.3 六個特徵的相關係數矩陣

以下用 python 實現邏輯回歸。實現了梯度下降、 IRLS 和通用牛頓法三種優化演算法。通用牛頓法是指採用 [3.7] 式對 $eta$ 進行更新。此時可以通過 $eta$ 控制學習速率。梯度下降法也用到了學習速率 $eta$ 。最後在測試集上查看模型的表現。此實現使用了 Numpy 的矩陣運算功能，使用了 Scikit-Learn 的數據標準化和模型評價功能。

# -*- coding: utf-8 -*-"""Created on Tue Feb 13 12:33:54 2018@author: zhangjuefei"""import numpy as npfrom sklearn.model_selection import train_test_splitfrom sklearn.preprocessing import StandardScalerfrom sklearn.metrics import classification_report, roc_curve, roc_auc_scoreimport matplotlib.pyplot as pltclass LogisticRegression: def __init__(self, method="IRLS", max_iters=10, eta=0.1): self.max_iters = int(max_iters) # max number of iterations self.method = str(method) # IRLS, Newton or Gradient self.eta = float(eta) # learning rate for Newton and Gradient method def fit(self, X, y): X = np.mat(X) y = np.mat(y).T self.beta = np.mat(np.random.random((X.shape[1], 1)) / 100.0) for i in range(self.max_iters): p = 1.0 / (1.0 + np.power(np.e, -np.matmul(X, self.beta))) p[p >= 1.0] = 1.0 - 1e-10 p[p <= 0.0] = 1e-10 # truncate the propabilities to prevent W from being singular W = np.mat(np.diag(np.array(np.multiply(p, 1.0-p)).ravel())) if self.method == "IRLS": z = X * self.beta + W.I * (y - p) self.beta = (X.T * W * X).I * X.T * W * z elif self.method == "Newton": self.beta = self.beta + self.eta * (X.T * W * X).I * X.T * (y - p) elif self.method == "Gradient": self.beta = self.beta + self.eta * X.T * (y - p) else: raise Exception("Unsupported method {}".format(self.method)) def predict(self, X): X = np.mat(X) p = 1.0 / (1.0 + np.power(np.e, -np.matmul(X, self.beta))) return (p > 0.5).astype(np.int), p X = np.loadtxt("D:/study/column_2C_weka.csv", skiprows=1, usecols=np.arange(0,6), delimiter=",")y = np.loadtxt("D:/study/column_2C_weka.csv", skiprows=1, usecols=6, delimiter=",", dtype=np.str)y = (y == "Abnormal").astype(np.int)train_x, test_x, train_y, test_y = train_test_split(X, y, train_size=0.6)ss = StandardScaler()train_x = ss.fit_transform(train_x)test_x = ss.transform(test_x)lr = LogisticRegression("Gradient", 1000, 0.01)lr.fit(train_x, train_y)predict, predict_proba = lr.predict(test_x)print(classification_report(test_y, predict))fpr, tpr, thresholds = roc_curve(test_y, predict_proba)fig, ax = plt.subplots(1, 1, figsize=(6,6))ax.plot(fpr, tpr)ax.plot([0.0, 1.0], [0.0, 1.0], "r--")ax.grid(True)ax.set_xlabel("False Positive Rate")ax.set_ylabel("True Positive Rate")ax.text(0.82, 0.1, "auc: {:.2f}".format(roc_auc_score(test_y, predict_proba)))for f, t, th in zip(fpr[::3], tpr[::3], thresholds[::3]): ax.text(f+0.01, t-0.03, "{:.3f}".format(th), rotation=-15, size=8)

模型表現：

precision recall f1-score support Normal 0.57 0.95 0.71 38 Abnormal 0.97 0.69 0.80 86avg / total 0.85 0.77 0.78 124

圖 4.4 模型 roc 曲線，曲線上的數字為概率閾值。

圖 4.5 六個特徵的權重

五、參考書籍

統計學習基礎(第2版)(英文) (豆瓣)book.douban.com

最優化導論 (豆瓣)book.douban.com