逻辑回归——牛顿法矩阵实现方式

引言

逻辑回归常用来处理分类问题，最常用来处理二分类问题。

生活中经常遇到具有两种结果的情况（冬天的北京会下雪，或者不会下雪；暗恋的对象也喜欢我，或者不喜欢我；今年的期末考试会挂科，或者不会挂科……）。对于这些二分类结果，我们通常会有一些输入变量，或者是连续性，或者是离散型。那么，我们怎样来对这些数据建立模型并且进行分析呢？

我们可以尝试构建一种规则来根据输入变量猜测二分输出变量，这在统计机器学上被称为分类。然而，简单的给一个回答“是”或者“不是”显得太过粗鲁，尤其是当我们没有完美的规则的时候。总之呢，我们不希望给出的结果就是武断的“是”或“否”，我们希望能有一个概率来表示我们的结果。

一个很好的想法就是，在给定输入\(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.