机器学习（7）——Logistic回归（从疝气病症预测病马的死亡率）

转载请注明作者和出处：https://blog.csdn.net/qq_28810395
Python版本： Python3.x
运行平台： Windows 10
IDE： Pycharm profession 2019
如有需要数据集或整个工程文件的可以关注Stefan0704回复Logistic回归

一、前言

我们日常生活中不能避免的要遇到很多最优化问题，比如以最小工作量却获得最大利益？如何在最短时间内从A点到B点？如何设计发动机使得油耗最小而功率最大?由此解决这些问题含有重要意义，也有现实需要。

二、Logistic回归

Logistic回归是一种广义线性回归，具体作用与意义不具体描述，想了解请看Logistic 回归详解，在这只是利用其进行分类，其主要思想为：更具现有数据对分类边界线建立回归公式，以此进行分类，Logistic回归可以是二分类的，也可以是多分类的，但是二分类的更为常用，也更加容易解释，多类可以使用softmax方法进行处理。实际中最为常用的就是二分类的Logistic回归。
分类与回归是机器学习的主要核心，回归的含义简单来说就是用一条直线（该线称为最佳拟合直线）对这些点进行拟合，这个拟合的过程就被称作回归。

优点：计算代价不高，易于理解和实现，存储资源低，适用数值型和标称型数据。
缺点：容易欠拟合，分类精度可能不高。

三、Logistic回归的一般过程

收集数据：采用任意方法收集数据。
准备数据：由于需要进行距离计算，因此要求数据类型为数值型。另外，结构化数据格式则最佳。
分析数据：采用任意方法对数据进行分析。
训练算法：大部分时间将用于训练，训练的目的是为了找最佳的分类回归系数。
测试算法：一旦训练步骤完成，分类将会很快。
使用算法：首先，我们需要输入一些数据，并将其转话换成对应的结构化数值；接着，基于训练好的回归系数就可以对这些数值进行简单的回归计算，判定他们属于那个类别；在这之后，我们就可以在输出的类别上做一些其他分析工作。

四、基于Logistic回归和Sigmoid函数的分类

对目标进行分类，就要求对负无穷到正无穷的输入，都可以进行分类，输出bool行数值（0或1），学习电信的或者自控的肯定接触过这种性质的函数，在这里不卖关子，就是我们学过的海维赛德阶跃函数，也叫做单位阶跃函数。由于海维赛德阶跃函数跳跃过程很难处理，为此在此介绍另一个阶跃函数——Sigmoid函数，数学公式为：
σ(z)=11+e−z\sigma (z) = \frac{1}{{1 + {e^{ - z}}}}σ(z)=1+e−z1
由图所示，输入随着xxx的逐渐增大，输出逐渐接近1，xxx逐渐减小，输出逐渐接近0，在xxx输入为0的时候，输出为0.5，分类将输出小于0.5的归为0类，大于0.5的归为1类。

确定了分类器函数形式后，就要确定最佳回归系数，下述进行讲解。

五、基于最优化方法的最佳回归系数确定

Sigmoid 函数的输入记为 z：由下面公式得出：
z=w0x0+w1x1+w2x2+⋅⋅⋅+wnxnz = {w_0}{x_0} + {w_1}{x_1} + {w_2}{x_2} + \cdot \cdot \cdot + {w_n}{x_n}z=w0x0+w1x1+w2x2+⋅⋅⋅+wnxn
采用向量写法：
z=wTxz = {w^T}xz=wTx
表示将这两个数值向量对应元素相乘然后全部加起来得到zzz值。向量xxx是分类器的输入数据，向量www也就是我们要找的最佳参数，从而使得分类器尽可能准确，

1.梯度上升法

梯度上升法基于的思想是：要找到某函数的最大值，最好的方法是沿着该函数的梯度方向探寻。如果梯度记为∇\nabla∇，则函数f(x,y)f(x,y)f(x,y)的梯度由下述公式表示：
▽f(x,y)=(∂f(x,y)∂x∂f(x,y)∂y)\triangledown f(x,y)=\binom{\frac{\partial f(x,y)}{\partial x}}{\frac{\partial f(x,y)}{\partial y}}▽f(x,y)=(∂y∂f(x,y)∂x∂f(x,y))

如上图所示，梯度上升算法到达一个点后重新估计移动的方向。P0−>P1,P1−>P2P0->P1,P1->P2P0−>P1,P1−>P2，根据梯度算子循环迭代，直到满足停止条件（迭代次数达到某个指定值或算法达到某个可以允许的误差范围）。迭代公式如下，ααα为步长，是移动量的大小。
w:=w+α▽f(w)w:=w+\alpha\triangledown{f(w)}w:=w+α▽f(w)
我们实际中经常用的为梯度向下，其公式如下：
w:=w−α▽f(w)w:=w-\alpha\triangledown{f(w)}w:=w−α▽f(w)

学习完理论，下面我们做个基于Logistic回归分类器的实验。

2.准备数据集

其TXT文件在公众号可获取。

3.训练算法：使用梯度上升找到最佳参数

在上述的数据集中，100个样本点中的每个样本点包含两个特征X1X1X1、X2X2X2,我们利用梯度上升法找其最佳回归系数，也就是拟合出Logistic回归模型的最佳参数。

梯度上升的伪代码如下：

每个回归系数初始化为1
重复R次
计算整个数据集的梯度
使用alpha×gradient更新回归系数的向量
返回回归系数

代码实现

import matplotlib.pyplot as plt
import numpy as np'''
函数说明：加载数据
Parameters:None
Returns:dataMat - 数据列表labelMat - 标签列表
'''
def loadDataSet():dataMat = []  # 创建数据列表labelMat = []  # 创建标签列表fr = open('testSet.txt')  # 打开文件for line in fr.readlines():lineArr = line.strip().split()  # 去回车，放入列表dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])  # 添加数据labelMat.append(int(lineArr[2]))  # 添加标签fr.close()  # 关闭文件return dataMat, labelMat
'''
函数说明:sigmodi函数
Paremeters:inX - 数据
Returns:sigmoid函数
'''
def sigmoid(inX):return 1.0/(1 + np.exp(-inX))'''
函数说明：梯度上升算法
Parameters:dataMatIn - 数据集classLables - 数据标签
Returns：weights.getA() - 求得的权重数组（最优参数）
'''
def gradAscent(dataMatIn, classLables):dataMatrix = np.mat(dataMatIn)  #转换成numpy的matlabelMat =  np.mat(classLables).transpose() #转换成numpy的mat，并进行转置m, n =np.shape(dataMatrix)#返回dataMatrix的大小。m为行，n为列alpha = 0.001  #移动补偿，也就是学习速率，控制更新的幅度maxCycles = 500 #最大迭代次数weights = np.ones((n,1))for k in range(maxCycles):h = sigmoid(dataMatrix *weights) #梯度上升矢量公式#权重系数计算公式error = labelMat - hweights = weights + alpha * dataMatrix.transpose()*errorreturn weights.getA()  #将矩阵转换为数组，返回权重数组

测试代码

import logRegres
import matplotlib.pyplot as plt
from numpy import *if __name__ == '__main__':dataArr, labelMat = logRegres.loadDataSet()print(logRegres.gradAscent(dataArr, labelMat))

结果

[[ 4.12414349][ 0.48007329][-0.6168482 ]]

4.分析数据：画出决策边界

我们已经解出了一组回归系数，他确定了不同类别数据之间的分割线，那么下面我就画出这条分割线，使优化过程更便于理解。

代码实现

 '''函数说明：绘制数据集和Logistic回归最佳拟合直线的函数Parameters:NoneReturns:None'''
def plotDataSet(weights):dataMat, labelMat = loadDataSet()  # 加载数据集dataArr = np.array(dataMat)  # 转换成numpy的array数组n = np.shape(dataMat)[0]  # 数据个数,即行数xcord1 = [] ; ycord1 = []  # 正样本xcord2 = [] ; ycord2 = []  # 负样本for i in range(n):if int(labelMat[i]) == 1: #1为正样本xcord1.append(dataMat[i][1])ycord1.append(dataMat[i][2])else:                     #0为负样本xcord2.append(dataMat[i][1])ycord2.append(dataMat[i][2])fig = plt.figure()ax = fig.add_subplot(111)   #添加subplotax.scatter(xcord1,ycord1,s=20,c='red',marker = 's', alpha=.5,label ='1') #绘制正样本ax.scatter(xcord2,ycord2,s=20,c='green',marker = 's', alpha=.5,label ='0') #绘制正样本x = np.arange(-3.0,3.0,0.1)y = (-weights[0] - weights[1] * x) / weights[2]ax.plot(x,y)plt.title('DataSet') #绘制titleplt.xlabel('x'); plt.ylabel('y') #绘制labelplt.legend()plt.show()

测试代码

if __name__ == '__main__':dataArr, labelMat = logRegres.loadDataSet()wigets = logRegres.gradAscent(dataArr, labelMat)logRegres.plotBestFit(wigets.getA())

绘图结果

其中绘制的分隔线设置了Sigmoid函数为0，回忆上一篇内容，0是两个分类的分界处，因此我们设定0=w0x0+w1x1+w2x2\large 0 = w_{0}x_{0}+w_{1}x_{1}+w_{2}x_{2}0=w0x0+w1x1+w2x2，然后接触X2X2X2和X1X1X1的关系式（及分割线的方程，注意X0=1X0=1X0=1）。
根据分类的结果来看，分类结果错分了几个点，效果已经很不错。但是尽管例子简单且数据量（100哥样本点）很小，却需要大量的计算（300次乘法），但真实情况下，数据集是非常巨大的，为此应对大数据集如果采用此算法计算成本太大。下面将对该算法进行改进，以计算更大的数据集与特征。

5.训练算法：随机梯度上升

由于梯度上升算法在每次更新回归系数时都需要遍历整个数据集，为此应对这个缺点，随机梯度上升算法的优势就体现出来了，算法思想是：一次仅用一个样本点来更新回归系数。

随机梯度上升算法伪代码：

所有回归系数初始化为1
对数据集中每个样本
计算该样本的梯度
使用alpha×gradient更新回归系数值
返回回归系数值

代码实现

#随机梯度上升法
def stocGradAscent0(dataMatrix, classLabels):m,n = shape(dataMatrix)alpha = 0.01weights = ones(n)   #initialize to all onesfor i in range(m):h = sigmoid(sum(dataMatrix[i]*weights))error = classLabels[i] - hweights = weights + alpha * error * dataMatrix[i]return weights

从程序代码看出，随机梯度上升算法与梯度上升算法很相似，但也存在区别：

后者的变量 h 和误差 error 都是向量，而前者则是数值；
前者没有矩阵的转换过程，所有变量的数据类型都是 numpy 数组。

测试代码

if __name__ == '__main__':dataArr, labelMat = logRegres.loadDataSet()wigets = logRegres.stocGradAscent0(array(dataArr), labelMat)logRegres.plotBestFit(wigets)

结果
由图看出，此次拟合的曲线错分了1/31/31/3的样本，效果相比于梯度上升算法欠佳，但这样直接比较有欠公平，一个算法时候优劣，核心比较是看它是否收敛，也就是说参数是否达到稳定值，是否还会存在不断变化。对此，我们在 stocGradAscent0() 函数上做了些修改，使其在整个数据集上运行 200 次。最终绘制的三个回归系数的变化情况。

由图看出，随机梯度上升算法在200迭代过程中，X2X2X2经过50次迭代后就达到稳定值，而X1X1X1、X0X0X0需要更多的迭代，另外值得注意，在迭代过程中还有一些小的周期波动，产生这种现象是由于存在一些不能正确分类的样本点（数据集非线性可分），为了弥补这个缺点，我们将对随机梯度上升算法进行改进。

6.算法改进：改进的随机梯度上升算法

代码实现
```
def stocGradAscent1(dataMatrix, classLabels, numIter=150):m,n = np.shape(dataMatrix) #返回dataMatrix的大小 m为行数 n为列数weights = np.ones(n) #[1. 1. 1.] 初始化参数 类型是array所以注意dataMatrix类型也应该是np.arrayfor j in range(numIter):dataIndex = list(range(m))  # [0,1,...99]for i in range(m):alpha = 4/(1.0+j+i)+0.01 #降低alpha的大小，每次减小1/(j+i)randIndex = int(random.uniform(0,len(dataIndex)))#随机选取样本h = sigmoid(sum(dataMatrix[randIndex]*weights))#选择随机选取的一个样本，计算h。对应位置相乘求和error = classLabels[randIndex] - hweights = weights+ alpha* error*dataMatrix[randIndex]#更新回归系数del(dataIndex[randIndex])#删除已经使用的样本return weights
```
该算法主要是在以下两方面做了改进。
- alpha在每次迭代的时候都会调整，随着迭代次数不断减小，但永远不会减小到0。这样做的原因是为了保证在多次迭代之后新数据仍然具有一定的影响。如果要处理的问题是动态变化的，那么可以适当加大上述常数项，来确保新的值获得更大的回归系数。而且在降低alpha的函数中，alpha每次减少1/(j+i)，其中j是迭代次数，i是样本的下标。这样当j<<max(i)时，alpha就不是严格下降的。避免参数的严格下降也常见于模拟退火算法等其他优化算法中。
- 这里通过每次随机选取一个样本来更新回归系数。这种方法可以有效的减少周期性的波动。

测试代码

if __name__ == '__main__':dataArr, labelMat = logRegres.loadDataSet()wigets  = logRegres.stocGradAscent1(array(dataArr), labelMat)logRegres.plotBestFit(wigets)

绘图结果
进行迭代比较

相比于梯度上升算法，随机梯度上升算法的收敛效果更好。我们一共有100个样本点，改进的随机梯度上升算法迭代次数为150.而上图显示150000次迭代次数的原因是，使用一次样本就更新一下回归系数，要迭代150次。所以150*100=150000次。迭代150次，更新1.5万次回归参数。从右图我们可以知道基本迭代2000次即迭代20次回归系数已经收敛了，比初始算法收敛要更快一些。

六、示例、从疝气病症预测病马的死亡率

收集数据：给定数据文件。
准备数据：用Python解析文本文件并填充缺失值。
分析数据：可视化并观察数据。
训练算法：使用优化算法，找到最佳的系数。
测试算法：为了量化回归的效果，需要观察错误率。根据错误率决定是否退到训练阶段，通过改变迭代的次数和步长等参数来得到更好的回归系数。
使用算法：实现一个简单的命令行程序来收集妈的症状并输出预测结果。

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

假设有100个样本和20个特征，这些数据都是机器手机回来的。若机器上的某个传感器损坏导致一个特征无效时该怎么办？此时是否要扔掉整个数据？这种情况下，另外19个特征怎么办？是否还可用？答案是肯定的。因为有的数据相当昂贵，扔掉和重新获取都是不可取的，下面给出了一些处理方法。

使用可用特征的均值来填补缺失值
使用特殊值来填补缺失值，如-1
忽略有缺失值的样本
使用相似样本的均值填补缺失值
使用另外的机器学习算法预测缺失值
现在我们要进行数据预处理。有两件事：
- 一、所有的缺失值必须用一个实数值；因为我们使用的Numpy数据类型不允许包含缺失值。这里选择实数0来替换所有缺失值。这样做的好处在于更新时不会影响回归系数的值，如果选择实数0来替换缺失值，那么dataMatrix的某特征对应值为0，那么该特征的系数将不做更新。所有以0代替缺失值在Logistic回归中可以满足这个要求。
- 二、如果在测试数据集中发现一条数据的类别标签已经缺失，那么我们的简单做法是将该条数据丢弃，这是因为类别标签与特征不同，很难确定采用某个合适的值来替换。采用Logistic回归进行分析时这类做法是合理的。现在我没有一个干净可用的数据集和一个不错的优化算法，将二者结合，预测马的生死问题。

2.测试代码：用Logistic回归进行分类

根据训练集的每个特征用于最优化算法中得到回归系数。把测试集的特征乘以回归系数并求和放入Sigmoid函数中。如果返回值大于0.5，预测标签为1，否则为0.

代码实现

def classifyVector(inX, weights):prob = sigmoid(sum(inX*weights))if prob > 0.5: return 1.0else: return 0.0def colicTest():frTrain = open('horseColicTraining.txt'); frTest = open('horseColicTest.txt')trainingSet = []; trainingLabels = []for line in frTrain.readlines():currLine = line.strip().split('\t')lineArr =[]for i in range(21):lineArr.append(float(currLine[i]))trainingSet.append(lineArr)trainingLabels.append(float(currLine[21]))trainWeights = stocGradAscent1(array(trainingSet), trainingLabels, 1000)errorCount = 0; numTestVec = 0.0for line in frTest.readlines():numTestVec += 1.0currLine = line.strip().split('\t')lineArr =[]for i in range(21):lineArr.append(float(currLine[i]))if int(classifyVector(array(lineArr), trainWeights))!= int(currLine[21]):errorCount += 1errorRate = (float(errorCount)/numTestVec)print ("the error rate of this test is: %f" % errorRate)return errorRatedef multiTest():numTests = 10; errorSum=0.0for k in range(numTests):errorSum += colicTest()print ("after %d iterations the average error rate is: %f" % (numTests, errorSum/float(numTests)))

测试代码

if __name__ == '__main__':logRegres.multiTest()

训练结果

the error rate of this test is: 0.328358
the error rate of this test is: 0.373134
the error rate of this test is: 0.343284
the error rate of this test is: 0.417910
the error rate of this test is: 0.328358
the error rate of this test is: 0.358209
the error rate of this test is: 0.313433
the error rate of this test is: 0.268657
the error rate of this test is: 0.388060
the error rate of this test is: 0.388060
after 10 iterations the average error rate is: 0.350746

由结果看来，在30%的数据缺失情况下，10次迭代之后的平均错误率是35%左右，已经很不错，修改迭代次数和步长可以更加降低错误率。

七、总结

Logistic回归的目的是寻找一个非线性函数Sigmoid的最佳拟合参数，求解过程可以由最优化算法来完成。
随机梯度上升算法与梯度上升算法的效果相当，但占用更少的计算资源。此外随机梯度算法是一个在线算法，可在新数据来时就完成参数更新，不需要重新读取整个数据集来进行批处理运算。

参考信息

https://www.jianshu.com/p/c9383ae24756

https://blog.csdn.net/qq_25174673/article/details/83749836

https://blog.csdn.net/qq_25174673/article/details/84069585