邏輯回歸——牛頓法矩陣實現方式

引言

邏輯回歸常用來處理分類問題，最常用來處理二分類問題。

生活中經常遇到具有兩種結果的情況（冬天的北京會下雪，或者不會下雪；暗戀的對象也喜歡我，或者不喜歡我；今年的期末考試會掛科，或者不會掛科……）。對於這些二分類結果，我們通常會有一些輸入變量，或者是連續性，或者是離散型。那么，我們怎樣來對這些數據建立模型並且進行分析呢？

我們可以嘗試構建一種規則來根據輸入變量猜測二分輸出變量，這在統計機器學上被稱為分類。然而，簡單的給一個回答“是”或者“不是”顯得太過粗魯，尤其是當我們沒有完美的規則的時候。總之呢，我們不希望給出的結果就是武斷的“是”或“否”，我們希望能有一個概率來表示我們的結果。

一個很好的想法就是，在給定輸入\(X\)的情況下，我們能夠知道Y的條件概率\(Pr(Y|X)\)。一旦給出了這個概率，我們就能夠知道我們預測結果的准確性。

讓我們把其中一個類稱為1，另一個類稱為0。（具體哪一個是1，哪一個是0都無所謂）。\(Y\)變成了一個指示變量，現在，你要讓自己相信，\(Pr(Y=1)=EY\)，類似的，\(Pr(Y=1|X=x)=E[Y|X=x]\)。

假設\(Y\)有10個觀測值，分別為 0 0 0 1 1 0 1 0 0 1.即6個0,4個1.那么，\(Pr(Y=1)=\frac{count(1)}{count(n)}=\frac{4}{10}=0.4\)，同時，\(EY=\frac{sum(Y)}{count(n)}=\frac{4}{10}=0.4\)

換句話說，條件概率是就是指示變量（即\(Y\))的條件期望。這對我們有幫助，因為從這個角度上，我們知道所有關於條件期望的估計。我們要做的最直接的事情是挑選出我們喜歡的平滑器，並估計指示變量的回歸函數，這就是條件概率函數的估計。

有兩個理由讓我們放棄陷入上述想法。

1. 概率必須介於0和1之間，但是我們在上面估計出來的平滑函數的輸出結果卻不能保證如此，即使我們的指示變量\(y_i\)不是0就是1；
2. 另一種情況是，我們可以更好地利用這個事實，即我們試圖通過更顯式地模擬概率來估計概率。

假設\(Pr(Y=1|X=x)=p(x;\theta)\),\(p\)是參數為\(\theta\)的函數。進一步，假設我們的所有觀測都是相互獨立的，那么條件似然函數可以寫成：

\[\prod _{i=1}^nPr(Y=y_i|X=x_i)=\prod _{i=1}^np(x_i;\theta)^{y_i}(1-p(x_i;\theta))^{1-y_i} \]

回憶一下，對於一系列的伯努利試驗\(y_1,y_2,\cdots,y_n\)，如果成功的概率都是常數\(p\)，那么似然概率為：

\[\prod _{i=1}^n p^{y_i}(1-p)^{1-y_i} \]

我們知道，當\(p=\hat{p}=\frac{1}{n}\sum _{i=1}^ny_i\)時，似然概率取得最大值。如果每一個試驗都有對應的成功概率\(p_i\)，那么似然概率就變成了

\[\prod _{i=1}^n p_i^{y_i}(1-p_i)^{1-y_i} \]

不添加任何約束的通過最大化似然函數來估計上述模型是沒有意義的。當\(\hat{p_i}=1\)的時候，\(y_i=1\)，當\(\hat{p_i}=0\)的時候，\(y_i=0\)。我們學不到任何東西。如果我們假設所有的\(p_i\)不是任意的數字，而是相互連接在一起的，這些約束給我們提供了一個很重要的參數，我們可以通過這個約束來求得似然函數的最大值。對於我們正在討論的這種模型，約束條件就是\(p_i=p(x_i;\theta)\)，當\(x_i\)相同的時候，\(p_i\)也必須相同。因為我們假設的\(p\)是未知的，因此似然函數是參數為\(\theta\)的函數，我們可以通過估計\(\theta\)來最大化似然函數。

邏輯回歸

總結一下：我們有一個二分輸出變量\(Y\)，我們希望構造一個關於\(x\)的函數來計算\(Y\)的條件概率\(Pr（{Y=1|X=x}）\)，所有的未知參數都可以通過最大似然法來估計。到目前為止，你不會驚訝於發現統計學家們通過問自己“我們如何使用線性回歸來解決這個問題”。

• 最簡單的一個想法就是令\(p(x)\)為線性函數。無論\(x\)在什么位置，\(x\)每增加一個單位，\(p\)的變化量是一樣的。由於線性函數不能保證預測結果位於0和1之間，因此從概念上線性函數就不適合。另外，在很多情況下，根據我們的經驗可知，當\(p\)很大的時候，對於\(p\)的相同的變化量，\(x\)的變化量將會大於\(p\)在1/2附近的變化量。線性函數做不到這樣。
• 另一個最直觀的想法是令\(log\ p(x)\)為\(x\)的線性函數。但是對數函數只在一個方向上是無界的，因此也不符合要求。
• 最后，我們選擇了\(p\)經過logit變換以后的函數\(ln\frac{p}{1-p}\)。這個函數就很好啊，滿足了我們的所有需求。

最后一個選擇就是我們所說的邏輯回歸。一般的，邏輯回歸的模型可以表示為如下形式：

\[ln\frac{p(x)}{1-p(x)}=\beta_0+x\beta \]

根據上式，解出\(p\)

\[p(x;b,w)=\frac{e^{\beta_0+x\beta}}{1+e^{\beta_0+x\beta}}=\frac{1}{1+e^{-(\beta_0+x\beta)}} \]

為了最小化錯分率，當\(p\ge 0.5\)的時候，我們預測\(Y=1\)，否則\(Y=0\)。這意味着當\(\beta_0+x\beta\)非負的時候，預測結果為1，否則為預測結果為0.因此，邏輯回歸為我么提供了一個線性分類器。決策邊界是\(\beta_0+x\beta=0\)，當\(x\)是一維的時候，決策邊界是一個點，當\(x\)是二維的時候，決策邊界是一條直線，以此類推。空間中某個點到決策邊界的距離為\(\beta_0/||\beta||+x\cdot\beta/||\beta||\).邏輯回歸不僅告訴我們兩個類的決策邊界，還以一種獨特的方式根據樣本點到決策邊界的距離給出該點分屬於某類的概率。當\(||\beta||\)越大的時候，概率取極端值（0或1）的速度就越快。上述說明使得邏輯回歸不僅僅是一個分類器，它能做出更簡健壯、更詳細的預測，並能以一種不同的方式進行擬合;但那些強有力的預測可能是錯誤的。

似然函數和邏輯回歸

因為邏輯回歸的預測結果是概率，而不是類別，因此我們可以用似然函數來擬合模型。對於每一個樣本點，我們有一個特征向量\(x_i\)，這個向量的維度就是特征的個數。同時還有一個觀測類別\(y_i\)。當\(y_i=1\)的時候，該類的概率為\(p\)，否則為\(1-p\)。因此，似然函數為：

\[L(\beta_0,\beta) = \prod _{i=1}^np(x_i)^{y_i}(1-p(x_i))^{1-y_i} \]

對數似然函數為：

\[\begin{aligned} \ell(\beta_0,\beta) &= \sum_{i=1}^ny_iln(p(x_i))+(1-y_i)ln(1-p(x_i)) \\&=\sum_{i=1}^n(y_i(ln\frac{p(x_i)}{1-p(x_i)})+ln(1-p(x_i))) \\&=\sum_{i=1}^ny_i(\beta_0+x_i\beta)-ln(1+e^{\beta_0+x_i\beta}) \end{aligned} \]

為了表示方便，統一將\(\beta_0,\beta\)表示成\(\beta\),則\(\ell\)對\(\beta\)的一階導數為：

\[\begin{aligned} \frac{\partial \ell}{\partial \beta} &=\sum _{i=1}^n[y_i-\frac{e^{\beta^Tx_i}}{1+\beta^Tx_i}]x_i\\& = \sum _{i=1}^n(y_i-p_i)x_i\end{aligned} \]

多分類邏輯回歸

如果\(Y\)有多個類別，我們仍然可以使用邏輯回歸。假如有\(k\)個類別，分別是\(0,1,\cdots,k-1\)，對於每一個類\(k\)，其都有對應的\(\beta_0\)和\(\beta\)，每個類對應的概率為:

\[P(Y=c|X=x)=\frac{e^{\beta_0^c+x\beta^c}}{\sum e^{\beta_0^c+x\beta^c}} \]

觀察上式可以發現，二分類邏輯回歸求是多分類邏輯回歸的特例.

在這里，讀者可能比較好奇，根據上式，二分類邏輯回歸的分母中的1是怎么來的。其實，無論有多少個類，我們總是將第一類的系數設置為0，那么類別為0的那部分在分母中對應的就是1.這樣做對模型的通用性沒有任何影響。

有讀者可能會問，為什么偏要把第一個類的系數設置為0，而不是其他的類。事實上，你可以設置任何一個類的系數為0，並且最終計算出來的結果都是一樣的。所以，按照慣例，我們都是把第一個類的系數設置為0.

牛頓法求解參數

為了求出待估參數\(\beta\)，我們利用Newton-Raphson算法。首先對對數似然函數求二階偏導：

\[\frac{\partial ^2\ell (\beta)}{\partial \beta \partial \beta^T}=-\sum_{i=1}^nx_ix_i^Tp_i(1-p_i) \]

注意，上面的\(x_i\)是個向量，也就是上面所說的特征向量，維度為特征個數加一。即假設原始數據為\(n\times m\)矩陣，其中n表示觀測數，m表示特征數。則\(x_i\)的長度為m+1。根據上述說明，上面的二階偏導實際上是一個\((m+1)\times (m+1)\)的矩陣。

如果給定一個\(\hat{\beta}^{old}\)，則一步牛頓迭代為（梯度下降）：

\[\hat{\beta}^{new}=\hat{\beta}^{old}-(\frac{\partial ^2\ell (\beta)}{\partial \beta \partial \beta^T})^{-1} \cdot\frac{\partial \ell({\beta})}{\partial \beta} \]

將上述式子表示成矩陣的形式就是：

\[\frac{\partial \ell(\beta)}{\partial \beta}=X^T(y-p) \]

\[\frac{\partial ^2\ell (\beta)}{\partial \beta \partial \beta^T}=-X^TWX \]

其中，\(X\)為原始自變量矩陣，\(y\)為類別向量，\(p\)為預測概率向量，\(W\)是一個\(n\times n\)對角矩陣，第\(i\)個元素取值為\(p(x_i,\hat{\beta}^{old})(1-p(x_i,\hat{\beta}^{old}))\).

聯立上述兩個式子，可以得出參數的迭代公式：

\[\begin{aligned}\hat{\beta}^{new} &=\hat{\beta}^{old}+(X^TWX)^{-1}X^T(y-p) \\ &=(X^TWX)^{-1}X^TW(X\hat{\beta^{old}}+W^{-1}(y-p)) \\ &=(X^TWX)^{-1}X^TWz \end{aligned} \]

其中，\(z=X\hat{\beta^{old}}+W^{-1}(y-p)\).

實際上，\(X^TWX\)是一個黑塞矩陣

\[H(\beta) = \frac{\partial ^2\ell (\beta)}{\partial \beta \partial \beta^T} \]

即目標函數對\(\beta\)的二階偏導，那么，上述迭代公式也可以寫作：

\[\begin{aligned} \hat{\beta}^{new} &=\hat{\beta}^{old} - H(\beta)^{-1}\frac{\partial \ell(\beta)}{\beta} \\ &= \hat{\beta}^{old} - H(\beta)^{-1}X^T(y-p) \end{aligned} \]

上述是我們熟悉的牛頓迭代公式。

矩陣相乘的計算不算復雜，但是當數據量上升以后，黑塞矩陣的求逆就非常復雜了，因此衍生出許多擬牛頓算法，本節不討論優化算法。

很明顯，本例的目標函數就是對數似然函數\(\ell(\beta)\)，也就是求其最大值。然而，很多同學已經習慣了牛頓法求最小值，因此，為了大家看着方便，下面介紹梯度下降法求解邏輯回歸。

只需要在上述似然函數前面加一個負號，本例就變成了一個梯度下降的問題了。為了形式上好看，還可以在前面對數似然函數求一個均值，即除以樣本量。

假設\(J(\beta)\)是我們的目標函數，則

\[\begin{aligned} J(\beta) &= -\frac{1}{n}\ell(\beta) \\ &=-\frac{1}{n}\sum_{i=1}^ny_iln(p(x_i))+(1-y_i)ln(1-p(x_i)) \end{aligned} \]

此時我們的梯度公式就變成了：

\[\frac{\partial J(\beta)}{\partial \beta}=-\frac{1}{n}\sum _{i=1}^n(y_i-p_i)x_i=\frac{1}{n}\sum _{i=1}^n(p_i-y_i)x_i \]

我們的二階偏導數就變成了

\[\frac{\partial ^2J(\beta)}{\partial \beta\partial\beta^T} =\frac{1}{n} \sum_{i=1}^nx_ix_i^Tp_i(1-p_i) \]

那么，此時為了求得我們的回歸系數，即求使得\(J(\beta)\)最小的系數。牛頓迭代公式就變成了：

\[\hat{\beta}^{new} = \hat{\beta}^{old}-H^{-1}\nabla=\hat{\beta}^{old}-\frac{1}{n}(X^T\cdot diag(p)\cdot diag(1-p) \cdot X)^{-1}\cdot \frac{1}{n}X^T(p-y) \]

按照上述思路，編程實現邏輯回歸求解是比較簡單的。

#迭代函數
from numpy import *
def sigmoid(x):
    return 1.0/(1.0+exp(-x))
def LogitReg(x,y,tol = 0.001,maxiter = 1000):
    samples,features = x.shape  #分別表示觀測樣本數量和特征數量
    features += 1
    
    #全部轉換為矩陣
    xdata = array(ones((samples,features)))
    xdata[:,0] = 1
    xdata[:,1:] = x
    xdata = mat(xdata)  #sample行，features列的輸入
    
    y = mat(y.reshape(samples,1))  #label，一個長度為samples的向量
    
    #首先初始化beta，令所有的系數為1,生成一個長度為features的列向量
    beta = mat(zeros((features,1)))
    
    iternum = 0 #迭代計數器
    
    #計算初始損失
    
    loss0 = float('inf')
    J = []
    while iternum < maxiter:
        try:
            p = sigmoid(xdata*beta) #計算似然概率
            nabla = 1.0/samples*xdata.T*(p-y)   #計算梯度
            H = 1.0/samples*xdata.T*diag(p.getA1())* diag((1-p).getA1())*xdata  #計算黑塞矩陣
            
            loss = 1.0/samples*sum(-y.getA1()*log(p.getA1())-(1-y).getA1()*log((1-p).getA1())) #計算損失
            J.append(loss)
            beta =beta -  H.I * nabla  #更新參數
            iternum += 1 #迭代器加一
            if loss0 - loss < tol:
                break
            loss0 = loss
        except:
            H = H + 0.0001
            #通常當黑塞矩陣奇異的時候，將矩陣加上一個很小的常數。
            break
        
    return beta,J

#預測函數
def predictLR(data,beta):
    data = array(data)
    if len(data.shape) == 1:
        length = len(data)
        newdata = tile(0,length+1)
        newdata[0] = 1
        newdata[1:] = data
        newdata = mat(newdata)
        pass
    else:
        shape = data.shape
        newdata = zeros((shape[0],shape[1]+1))
        newdata[:,0] = 1
        newdata[:,1:] = data
        newdata = mat(newdata)
    return sigmoid(newdata*beta)

df = pd.read_csv('df.csv',header=None)
df = array(df)
df.shape
xdata = df[:,:3]
ydata = df[:,3]

beta,J = LogitReg(xdata,ydata) #擬合

testdata = xdata[1:10,]

predictLR(testdata,beta)

matrix([[ 0.4959212 ],
        [ 0.44642627],
        [ 0.47419207],
        [ 0.42209742],
        [ 0.41802565],
        [ 0.51283217],
        [ 0.44833226],
        [ 0.41252982],
        [ 0.47853786]])

http://www.stat.cmu.edu/~cshalizi/uADA/12/lectures/ch12.pdf

[1]周志華.機器學習[M].清華大學出版社,2016.

[2]李航著.統計學習方法[M].清華大學出版社,2012.