逻辑回归从理论到 Python 实现

1. 逻辑回归的模型函数

前面我们讲了线性回归模型，将线性模型用于回归问题中。这篇我们讲一下线性模型用于分类任务。在二分类问题中，对于线性回归所产生的预测值

z = \theta^Tx

我们需将这个预测值 $z$ 转化为0/1值，最理想的是“单位阶跃函数”，即若预测值 $z$ 大于零就判为1，若预测值 $z$ 小于零则判为反例，预测值为临界值则可以任意判别。但是由于单位阶跃函数不连续，我们希望可以找到一个连续的，同时在一定程度上近似“单位阶跃函数“的函数，那么逻辑斯特函数(logistic function)则是一个很不错的替代函数：

y = \frac{1}{1+e^{-z}}

如下图所示，这是一种"Sigmoid"函数，它将 $z$ 值转化为一个接近0或1的y值。

我们将 $z = \theta^Tx$ 代入上式，最终得到逻辑回归的模型函数，如下

y = \frac{1}{1+e^{-\theta^Tx}}

作者个人理解：“单位阶跃函数”是将预测值直接转换为分类，而“逻辑斯特函数”是将预测值转换为0-1之间的一个数y，这个值正好可以当成某一类别的概率值，从而实现分类。

2. 逻辑回归的目标函数

逻辑回归的模型函数已经知道，求逻辑回归的模型函数，其实就是确定该模型函数中 $\theta$ 这个参数的值。在讲如何求这个参数之前，需要求得逻辑回归的目标函数，而在求解目标函数的过程中，需先讲解一下以下知识点：

一个事件的几率（odds）：指该事件发生与不发生的概率比值，若事件发生概率为 $p$ ，那么事件发生的几率就是

odds = \frac{p}{1-p}

那么该事件的**对数几率（log odds,亦称logit）**就是：

logit(p) =ln \frac{p}{1-p}

将式子

y = \frac{1}{1+e^{-\theta^Tx}}

转化为如下

ln\frac{y}{1-y} = \theta^Tx

也就是说，输出y的对数几率是由输入x的线性函数表示的模型，这就是逻辑回归模型。当 $\theta^Tx$ 的值越接近正无穷，y的值是越接近1的，那么我们将y视为样本x为正例，即y=1的概率。那么上述公式可转化为如下

ln{\frac{p(y=1|x)}{p(y=0|x)}}=\theta^Tx

从而得

p(y=1|x)=\frac{e^{\theta^Tx}}{1+e^{\theta^Tx}}\\ p(y=0|x)=\frac{1}{1+e^{\theta^Tx}}

进一步可转化为

p(y|x) = y\frac{e^{\theta^Tx}}{1+e^{\theta^Tx}}+(1-y)\frac{1}{1+e^{\theta^Tx}}

在统计学中，常常使用极大似然估计法来求解，即找到一组参数，使得在这组参数下，我们的数据的似然度（概率）最大，似然函数如下：

L(\theta)=\Pi_{i=1}^{m} (y_i\frac{e^{\theta^Tx_i}}{1+e^{\theta^Tx_i}}+(1-y_i)\frac{1}{1+e^{\theta^Tx_i}})

取对数似然函数得

\begin{aligned} lnL(\theta)&=\sum_{i=1}^{m}ln(y_i\frac{e^{\theta^Tx_i}}{1+e^{\theta^Tx_i}}+(1-y_i)\frac{1}{1+e^{\theta^Tx_i}})\\ &=\sum_{i=1}^{m}(y_i\theta^Tx_i-ln(1+e^{\theta^Tx_i})) \end{aligned}

我们要最大化上式，也就等价于最小化

J(\theta)=-lnL(\theta)=\sum_{i=1}^{m}(-y_i \theta^Tx_i+ln(1+e^{\theta^Tx_i}))

上式即可当做逻辑回归的目标函数

在逻辑回归模型中，我们最大化似然函数和最小化对数似然损失函数实际上是等价的。

3. 求解目标函数

上一块知识讲解中已经得到了目标函数，这一块知识讲解中，主要是讲解如何得到目标函数的最小值。自然而来的还是想到了梯度下降算法。我们求上述目标函数的梯度得

\nabla J(\theta) =(\frac{\partial J(\theta)}{\partial \theta_1},\frac{\partial J(\theta)}{\partial \theta_2},...,\frac{\partial J(\theta)}{\partial \theta_{n+1}})

其中

\begin{aligned} \frac{\partial J(\theta)}{\partial \theta_j}&=\sum_{i=1}^{m}(\frac{1}{1+e^{-\theta^Tx^{(i)}}}-y^{(i)})x_j^{(i)} \end{aligned}

那么对于 $\theta_j$ 的一次迭代如下:

\begin{aligned} \theta_j &= \theta_j-\alpha\frac{\partial J(\theta)}{\partial \theta_j}\\ &=\theta_j-\alpha\sum_{i=1}^{m}(\frac{1}{1+e^{-\theta^Tx^{(i)}}}-y^{(i)})x_j^{(i)}\\ \end{aligned}

从而推出对于 $\theta$ 的一次迭代如下（并进行了向量化）

\begin{aligned} \theta &= \theta-\alpha \sum_{i=1}^m(\frac{1}{1+e^{-\theta^Tx^{(i)}}}-y^{(i)})(x^{(i)})^T\\ &=\theta-\alpha X^TE \end{aligned}

其中

E =\begin{Bmatrix} \frac{1}{1+e^{-\theta^Tx^{(1)}}}-y^{(1)}\\ \frac{1}{1+e^{-\theta^Tx^{(2)}}}-y^{(2)}\\ ...\\ \frac{1}{1+e^{-\theta^Tx^{(m)}}}-y^{(m)} \end{Bmatrix}

上面将求和公式转化为了向量的形式，具体推导的过程已省略，可自行推导，也可以参考文末的梯度下降算法的向量化推导的参考链接。

4. 代码的实现

上面讲了逻辑回归的理论，下面我们通过代码来实现逻辑回归，例子主要是机器学习实战书上的例子。首先来看一下部分数据集（数据集的地址见文末给的链接地址），如下

-0.017612	14.053064	0
-1.395634	4.662541	1
-0.752157	6.538620	0
-1.322371	7.152853	0
...
-2.168791	0.143632	1
1.388610	9.341997	0
0.317029	14.739025	0

我们需要从文件中读取数据，并返回读取之后的列表结果，实现的代码如下

'''
从文件中读取数据
'''
def loadDataSet(filePath):
    dataList = []
    labelList = []
    f = open(filePath)
    for line in f.readlines():
        lineArr = line.strip().split()
        dataList.append([1.0, float(lineArr[0]), float(lineArr[1])])
        labelList.append(int(lineArr[2]))
    return dataList, labelList

之后我们采用梯度下降算法进行求解，在迭代过程中使用向量化，实现的代码如下：

'''
sigmoid函数
'''
def sigmoid(inX):
    return 1.0 / (1 + np.exp(-inX))


'''
梯度下降函数
'''
def gradAscent(dataList, labelList):
    dataMat = np.mat(dataList)  # x数据转化为矩阵
    labelMat = np.mat(labelList).T  # y数据转换为矩阵并转置为列向量
    m, n = np.shape(dataMat)  # 返回矩阵的大小
    alpha = 0.001  # 步长
    maxCycles = 500  # 迭代次数
    weights = np.ones((n, 1))  # 权重列向量
    # 梯度下降算法
    for k in range(maxCycles):
        h = sigmoid(dataMat * weights)
        error = h - labelMat
        weights = weights - alpha * dataMat.T * error
    return weights

通过函数gradAscent可以得到一组回归系数，它确定了不同类别数据之间的分隔线，下面我们将这个分割线画出来，代码如下：

'''
画分类的分割线
'''
def plotBestFit(dataList, labelList, weights):
    dataArr = np.array(dataList)
    m = np.shape(dataArr)[0]  # 获取样本的数量
    xcord0 = []
    ycord0 = []
    xcord1 = []
    ycord1 = []
    # 对不同标记的数据进行分类
    for i in range(m):
        if int(labelList[i]) == 0:
            xcord0.append(dataArr[i, 1])
            ycord0.append(dataArr[i, 2])
        else:
            xcord1.append(dataArr[i, 1])
            ycord1.append(dataArr[i, 2])
    # 画散点图
    plt.scatter(xcord0, ycord0, s=30, c='red', marker='s')
    plt.scatter(xcord1, ycord1, s=30, c='green')
    # 画分割线
    x1 = np.arange(-3.0, 3.0, 0.1)
    x2 = (-weights[0] - weights[1] * x1) / weights[2]
    plt.plot(x1, x2)
    plt.xlabel('X1')
    plt.ylabel('X2')
    plt.show()

上述代码中，比较令人不解的是x2 = (-weights[0] - weights[1] * x1) / weights[2]这个是怎么来的？因为对于样本数据来说

0=w_0x_0+w_1x_1+w_2x_2

是两个分类（类别1和类别0）的分界处，根据该式代入 $x_1$ 即可求 $x_2$ 从而这条线可以在图中表示出来。运行如下代码：

1
2
3

dataList, labelList = loadDataSet("testSet.txt") 
weights = gradAscent(dataList, labelList)  
plotBestFit(dataList, labelList, weights.getA())

输出结果如下图所示，从图中可以看到分类效果相当不错，从图上看值只分错了几个点。

然而跟线性回归中使用梯度下降算法求最小值类似，这种方法需要大量的计算，每次更新回归系数都需要遍历整个数据集，假如数据量一多，该方法的计算复杂度就很高了。所以我们对上述的批量梯度下降算法进行改变，改进方法是一次仅用一个样本点来更新回归系数，该方法称为随机梯度下降算法。下面是随机梯度下降算法的实现代码

'''
随机梯度下降算法
'''
def stocGradAscent0(dataList, labelList):
    dataArr = np.array(dataList)  # x数据转化为矩阵
    m, n = np.shape(dataArr)
    alpha = 0.01
    weights = np.ones(n)
    for i in range(m):
        h = sigmoid(np.sum(dataArr[i]*weights))
        error = (h-labelList[i])
        weights = weights-alpha*error*dataArr[i]
    return weights

运行如下代码，得到了该方法的分类效果图：

1
2
3

dataList, labelList = loadDataSet("testSet.txt")
weights = stocGradAscent0(dataList, labelList)
plotBestFit(dataList, labelList, weights)

相比批量梯度下降算法而言，效果较差。直接比较效果是有点不公平的，因为两者使用的数据集次数不同。而**判断一个优化算法优劣的可靠方法是看它是否收敛，也就是说参数是否达到稳定值。**通过下面这段代码，来绘制三个回归系数的变化情况：

'''
记录随机梯度下降过程中的历史权重值
'''
def stocGradAscent0WeightHistory(dataList, labelList):
    dataMat = np.array(dataList)  # x数据转化为矩阵
    m, n = np.shape(dataMat)
    alpha = 0.01
    times = 100
    weights = np.ones(n)
    weightsHistory = np.zeros((times*m, n) )
    for j  in range(times):
        for i in range(m):
            h = sigmoid(np.sum(dataMat[i] * weights))
            error = (h - labelList[i])
            weights = weights - alpha * error * dataMat[i]
            weightsHistory[j*m+i, :] = weights
    return weightsHistory
'''
绘制历史权值图
'''
def plotSDGError(weightsHistory):
    plt.rcParams['font.sans-serif'] = ['SimHei']  # 用来正常显示中文标签
    plt.rcParams['axes.unicode_minus'] = False  # 用来正常显示负号
    fig = plt.figure()
    ax = fig.add_subplot(311)
    ax.plot(weightsHistory[:,0])
    plt.ylabel('X0')
    ax = fig.add_subplot(312)
    ax.plot(weightsHistory[:,1])
    plt.ylabel('X1')
    ax = fig.add_subplot(313)
    ax.plot(weightsHistory[:,2])
    plt.ylabel('X2')
    plt.xlabel(u'迭代次数')
    plt.show()
    
    
def main():
    dataList, labelList = loadDataSet("testSet.txt")
    weightsHistory = stocGradAscent0WeightHistory(dataList, labelList)
    plotSDGError(weightsHistory)

if __name__ == '__main__':
    main()

最终绘制的三个回归系数的变化情况如下：

从图中可以看到在大的波动之后,有些系数还会存在一些小的周期性波动。这种现象的原因是因为存在一些不能正确分类的样本点，在每次迭代时会引发系数的剧烈改变。我们希望算法能避免来回波动，从而收敛于某个值，同时是加快收敛速度，因此我们将上述的随机梯度下降算法进行更改。更改为如下：

'''
改进的随机梯度下降算法
'''
def stocGradAscent1(dataList, labelList, numIter=150):
    dataArr = np.array(dataList)
    m,n = np.shape(dataArr)
    weights = np.ones(n)
    for j in range(numIter):
        dataIndex = list(range(m))
        for i in range(m):
            alpha = 4/(1.0+j+i)+0.01        # 步长为动态变化
            rand = int(np.random.uniform(0, len(dataIndex)))
            choseIndex = dataIndex[rand]
            h = sigmoid(np.sum(dataArr[choseIndex]*weights))
            error = h-labelList[choseIndex]
            weights = weights-alpha*error*dataArr[choseIndex]
            del(dataIndex[rand])
    return weights

相对第一个的随机梯度下降算法，上述算法一方面使alpha在每次迭代的时候都会进行调整，缓解了数据波动或者高频波动。虽然alpha会随着迭代次数不断减小，但是不会减少至0，而这样子做也保证了在多次迭代之后新数据仍然具有一定的影响。另一方面，通过随机选取样本来更新回归系数，从而减少周期性的波动。同时每次随机从列表中选出一个值之后，就从列表中删除该值。

上述算法对机器学习实战中的使用的算法有进行更改，作者个人认为书上在选出一个值并删除该值之后的实现方式不对。按照书上实现的算法，比如我们dataIndex序列是[0, 1, 2, 3, 4, 5]，那么假如第一次随机产生的数是0，那么使用的则是dataArr[0]这个样本，之后使用del删除了0下标这个数。那么dataIndex序列变成了[1, 2, 3, 4, 5]。之后再随机产生一个数，比如还是0，那么使用的还是dataArr[0]这个样本，之后del删除了下标为0的这个数，那么序列就变成了[2, 3, 4, 5]了。那么这种实现方式跟书的作者意思就有点不一致了。书中作者的意思应该是从[0, 1, 2, 3, 4, 5]这个列表中随机选取一个值，比如选取了0，然后使用第0个样本，之后删除0这个数，那么下次就变成从[1, 2, 3, 4, 5]这个列表中随机选取一个值，比如选取了1，然后使用第1个样本，之后删掉1这个数。

运行如下代码，得到的效果如下图所示，该效果达到了跟批量梯度下降算法差不多的效果，但是所使用的计算机更少

1
2
3

dataList, labelList = loadDataSet("testSet.txt")
weights = stocGradAscent1(dataList, labelList)
plotBestFit(dataList, labelList, weights)

使用更改之后的算法，来显示每次迭代时各个回归系数的变化情况，实现的代码如下

'''
记录改进的随机梯度下降过程中的历史权重值
'''
def stocGradAscent1WeightHistory(dataList, labelList, numIter=150):
    dataArr = np.array(dataList)
    m, n = np.shape(dataArr)
    weights = np.ones(n)
    weightsHistory = np.zeros((numIter*m, n))
    for j in range(numIter):
        dataIndex = list(range(m))
        for i in range(m):
            alpha = 4 / (1.0 + j + i) + 0.01  # 步长为动态变化
            rand = int(np.random.uniform(0, len(dataIndex)))
            choseIndex = dataIndex[rand]
            h = sigmoid(np.sum(dataArr[choseIndex] * weights))
            error = h - labelList[choseIndex]
            weights = weights - alpha * error * dataArr[choseIndex]
            weightsHistory[j*m+i, :] = weights
            del (dataIndex[rand])
    return weightsHistory

运行如下代码，得到的效果图如图所示。从图中可以看到没有出现周期性的波动，这归功于样本随机选择机制，同时收敛的速度也更快。

1
2
3

dataList, labelList = loadDataSet("testSet.txt")
weightsHistory = stocGradAscent1WeightHistory(dataList, labelList, 40)
plotSDGError(weightsHistory)

5. 实战：从疝气病症预测病马的死亡率

下面使用Logistic回归来预测患有疝病的马的存活问题，数据集中包含了368个样本和28个特征。数据集中包含了医院检测马疝病的一些指标，有的指标比较主观，有的指标难以测量，例如马的疼痛级别。另外需要说明的是，除了部分指标主观和难以测量外，该数据还存在一个问题，数据集中有30%的值是缺失的。首先在使用Logistic回归预测病马的生死之前，需要处理数据集中的数据缺失问题。

5.1. 准备数据：处理数据中的缺失值

对于有些存在缺失的数据来说，扔掉和重新获取是不可取的，所以有以下这些方法来解决数据缺失的问题：

使用可用特征的均值来填补缺失值
使用特殊值来填补缺失值，如-1
忽略有缺失值的样本
使用相似样本的均值添补缺失值
使用另外的机器学习算法预测缺失值

对于该实战中使用的数据集，在预处理阶段需要做两件事：

所有的缺失值必须用一个实数值来替换，这里选择实数0来替换所有缺失值，恰好能适用于Logistic回归，这样做在更新时不会影响回归系数的值。另外由于sigmoid(0)=0.5，即它对结果的预测不具有任何倾向性，因此上述做法也不会对误差造成任何影响。
测试数据集中发现一条数据的类别标签已经缺失，那么应将这条数据丢弃，这是因为类别标签与特征不同，很难确定采用某个合适的值来替换

机器学习中如何处理缺失数据这个问题没有标准答案，取决于实际应用中的需求。

原始的数据集经过预处理之后保存为了两个文件：horseColicTest.txt和horseColicTraining.txt。这两个数据集和原始数据集见文末给出的链接。

5.2 用Logistic回归进行分类

使用Logistic回归方法进行分类，所需要做的就是把测试集上的每个特征向量乘以最优化方法得来的回归系数，再将该乘积结果求和，最后输入到Sigmoid函数中，如果对应的函数值大于0.5就预测类别标签为1。下面看看实战中的代码，首先是文件内容的读取、Sigmoid函数以及梯度下降算法函数，这边不再赘述：

'''
数据加载
'''
def loadDataSet(filePath):
    f = open(filePath)
    dataList = []
    labelList = []
    for line in f.readlines():
        currLine = line.strip().split('\t')
        lineArr = []
        for i in range(21):
            lineArr.append(float(currLine[i]))
        dataList.append(lineArr)
        labelList.append(float(currLine[21]))
    return dataList, labelList


'''
sigmoid函数
'''
def sigmoid(inX):
    return 1.0 / (1 + np.exp(-inX))


'''
随机梯度下降算法
'''
def stocGradAscent(dataList, labelList, numIter=150):
    dataArr = np.array(dataList)
    m, n = np.shape(dataArr)
    weights = np.ones(n)
    for j in range(numIter):
        dataIndex = list(range(m))
        for i in range(m):
            alpha = 4 / (1.0 + j + i) + 0.01  # 步长为动态变化
            rand = int(np.random.uniform(0, len(dataIndex)))
            choseIndex = dataIndex[rand]
            h = sigmoid(np.sum(dataArr[choseIndex] * weights))
            error = h - labelList[choseIndex]
            weights = weights - alpha * error * dataArr[choseIndex]
            del (dataIndex[rand])
    return weights

接下去是分类函数，该函数以回归系数和特征向量作为输入来计算对应的Sigmoid值，如果Sigmoid值大于0.5则函数返回1。

'''
进行分类
'''
def classifyVector(inX, weights):
    prob = sigmoid(np.sum(inX * weights))
    if prob > 0.5:
        return 1.0
    else:
        return 0.0

再接下去是在测试集上计算分类错误率，这边返回的错误率多次运行得到的结果可能稍有不同，主要是因为其中有随机的成分在里面。只有当梯度下降算法得到的回归系数已经完全收敛，那么结果才是确定的。

'''
在测试集上计算分类错误率
'''
def colicTest(trainWeights, testDataList, testLabelList):
    errorCount = 0  # 判断错误的数量
    testCount = 0  # 测试集总共的数量
    testCount = len(testDataList)
    for i in range(testCount):
        if int(classifyVector(np.array(testDataList[i]), trainWeights)) != \
                int(testLabelList[i]):
            errorCount += 1
    errorRate = float(errorCount)/testCount
    return errorRate

最后只要调用如下的主函数即可获得10次运行之后的平均分类错误率

def main():
    numTests = 10
    errorSum = 0.0
    trainDataList, trainLabelList = loadDataSet("horseColicTraining.txt")
    testDataList, testLabelList = loadDataSet("horseColicTest.txt")
    for i in range(numTests):
        trainWeights = stocGradAscent(trainDataList, trainLabelList, 500)
        errorSum += colicTest(trainWeights, testDataList, testLabelList)
    print("%d次的平均错误率为%f"%(numTests, errorSum/numTests))

某次运行程序之后，得到的平均分类错误率如下（即输出结果），如果调整stocGradAscent函数中参数numIter的值和梯度下降算法中的步长，那么平均可以降到20%左右。

1	`10次的平均错误率为0.329851`

附：数据集及代码

本文所有数据集及代码，均可在https://github.com/DawnGuoDev/MachineLearningStudy/tree/master/LogisticRegression中获取。

本文参考与作者推荐

机器学习，周志华
机器学习实战，[美]Peter Harrington
CSDN_梯度下降法的向量化推导
CSDN_梯度上升算法与随机梯度上升算法的实现