目录
5.1 基于logistic回归和sigmoid函数的分类
5.2基于最优化方法的最佳回归系数确定
5.2.1 梯度上升法
5.2.2 训练算法：使用梯度上升法找到最佳参数
5.2.3 分析数据：画出决策边界
5.2.4 训练算法：随机梯度上升
5.3 示例：从疝气病症预测病马的死亡率
5.3.1 准备数据：处理数据中的缺失值
5.3.2 测试算法：用logistic回归进行分类
5.4 本章小结
假设现在有一些数据点，我们用一条直线对这些点进行拟合（改线称为最佳拟合直线），这个拟合过程就叫作回归。利用logistic回归进行分类的主要思想是：根据现有数据对分类边界线建立回归公式，以此进行分类。这里的“回归”一词源于最佳拟合，表示要找到最佳拟合参数集。训练分类器时的做法就是寻找最佳拟合参数，使用的是最优化算法。
logistic回归的一般过程
（1）收集数据：采用任意方法收集数据。
（2）准备数据：由于需要进行距离计算，因此要求数据类型为数值型。另外，结构化数据格式则最佳。
（3）分析数据：采用任意方式对数据进行分析。
（4）训练算法：大部分时间将用于训练，训练的目的是为了找到最佳的分类回归系数。
（5）测试算法：一旦训练步骤完成，分类将会很快。
（6）使用算法：首先，需要输入一些数据，并将其转换成对应的结构化数值；接着，基于训练好的回归系数就可以对这些数值进行简单的回归计算，判定它们属于哪个类别；在这之后，就可以在输出的类别上做一些其他分析工作。
5.1 基于logistic回归和sigmoid函数的分类
logistic回归
优点：计算代价不高，易于理解和实现。
缺点：容易欠拟合，分类精度可能不高。
适用数据类型：数值型和标称型数据。
我们想要的函数应该是，能接受所有的输入然后预测出类别。例如，在两个类的情况下，上述函数输出0或1，该函数称为‘海维赛德阶跃函数（heaviside step function）’或直接称为单位阶跃函数。然而，海维赛德阶跃函数的问题在于：该函数在阶跃点上从0瞬间跳跃到1，这个瞬间跳跃过程有时很难处理。幸好，另一个函数也有类似的性质，且数学上更易处理，这就是sigmoid函数。sigmoid函数具体的计算公式如下：
下图给出了sigmoid函数在不同坐标尺度下的两条曲线图。当x为0时，sigmoid函数值为0.5.随着x的增大，对应的sigmoid值将逼近于1；而随着x的减小，sigmoid值将逼近于0.如果横坐标刻度足够大，sigmoid函数看起来很像一个阶跃函数：
两种坐标尺度下的sigmoid函数图。上图的横坐标是-6到6，这时的曲线变化较为平滑；下图横坐标尺度足够大，可以看到，在x=0点处sigmoid函数看起来很像阶跃函数。
因此，为了实现logistic回归分类器，可以在每个特征上都乘于一个回归系数，然后把所有的结果值相加，将这个总和带入sigmoid函数中，进而得到一个范围在0~1之间的数值。任何大于0.5的数据被分为1类，小于0.5即被归入0类。所以，logistic回归也可以被看成是一种概率估计。
确定了分类器的函数形式之后，现在的问题变成了：最佳回归系数是多少？如何确定它们的大小？
5.2基于最优化方法的最佳回归系数确定
sigmoid函数的输入记为z，由下面公式得出：
如果采用向量的写法，上述公式可以写成
，它表示将这两个数值向量对应元素相乘然后全部加起来即得到z值。其中的向量x是分类器的输入数据，向量w就是要找的最佳参数（系数），从而使得分类器尽可能地精确。为了寻找该最佳参数，需要用到最优化理论的一些知识。
5.2.1 梯度上升法
梯度上升法基于的思想是：要找到某函数的最大值，最好的方法就是沿着该函数的梯度方向探寻。如果梯度记为
，则函数f（x,y）的梯度由下式表示：
这是机器学习中最易造成混淆的一个地方，但在数学上并不难，需要做的只是牢记这些符号的意义。这个梯度意味着要沿x的方向移动
,沿y的方向移动
。其中，函数f(x,y)必须要在待计算的点上有定义并且可微。
梯度上升算法到达每个点后都会重新估计移动的方向。从p0开始，计算完该点的梯度，函数就根据梯度移动到下一点p1.在p1点，梯度再次被重新计算，并沿新的梯度方向移动到p2.如此循环迭代，知道满足停止条件。迭代的过程中，梯度算子总是保证我们能选取到最佳的移动方向。
从上图中可以看到，梯度算子总是指向函数值增长最快的方向。这里所说的是移动方向，而未提到移动量的大小。该量值称为步长，记作
。用向量来表示的话，梯度算法的迭代公式如下：
该公式将一直被迭代执行，直至达到某个停止条件为止，比如迭代次数达到某个指定值或算法达到某个可以允许的误差范围。
梯度下降法：梯度下降法和梯度上升法是一样的。梯度上升法用来求函数的最大值，而梯度下降法用来求函数的最小值。
5.2.2 训练算法：使用梯度上升法找到最佳参数
上图中有100个样本点，每个点包含两个数值型特征：x1和x2.在此数据集上，我们将通过使用梯度上升法找到最佳回归系数，也就是拟合出logistic回归模型的最佳参数。
梯度上升法的伪代码：
每个回归系数初始化为1
重复R次：
。。计算整个数据集的梯度
。。使用alphaxgradient更新回归系数的向量
。。返回回归系数
下面用代码实现梯度上升法：
首先做一个包含100个样本点的数据集放到testSet.txt文件中（其中第一个数是x1，第二个数是x2，第三个数代表类别）：
创建名为logRegres.py的文件，输入下列代码：
from numpy import *
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]))
return dataMat,labelMat
def sigmoid(inX):
return 1.0/(1+exp(-inX))
def gradAscent(dataMatIn,classLabels):
dataMatrix=mat(dataMatIn)
labelMat=mat(classLabels).transpose()
m,n=shape(dataMatrix)
alpha=0.001
maxCycles=500
weights=ones((n,1))
for k in range(maxCycles):
h=sigmoid(dataMatrix*weights)
error=labelMat-h
weights=weights+alpha*dataMatrix.transpose()*error
return weights
程序开头提供了一个便利函数loadDataSet()，它的主要功能是打开文本文件testSet.txt并逐行读取。每行前两个值分别是x1和x2，第三个值是数据对应的类别标签。此外，为了方便计算，该函数还将x0值设为1.0。
梯度上升法的实际工作是在函数gradAscent()里完成的，该函数有两个参数。第一个参数是dataMatIn，它是一个2维的numpy数组，每列分别代表每个不同的特征，每行则代表每个训练样本。现采用100个样本的简单数据集，它包含了两个特征x1和x2，再加上第0维特征x0，所以dataMatIn里存放的将是100x3的矩阵。mat（）将获得的输入数据转换成numpy矩阵。第二个参数是类别标签，它是一个1x100的行向量。为了便于矩阵运算，需要将该行向量转换成列向量，做法是将原向量转置，再将它赋值给labelMat。
变量alpha是向目标移动的步长，maxCycle是迭代次数。在for循环迭代完成后，将返回训练好的回归系数。需要强调的是，h=sigmoid(dataMatrix*weights)中的运算是矩阵运算。变量h不是一个数而是一个列向量，列向量的元素个数等于样本个数，这里是100.对应地，运算dataMatrix*weights代表的不止一次乘积运算，事实上该运算包含了300次的乘积。
运行代码：
if __name__=="__main__":
dataArr,labelMat=loadDataSet()
print(gradAscent(dataArr,labelMat))
运行结果：
[[ 3.8218897 ]
[ 0.06112853]
[-0.5704236 ]]
5.2.3 分析数据：画出决策边界
上面解出的一组回归系数，确定了不同类别数据之间的分隔线，现画出分隔线，打开logRegres.py添加以下代码：
#画出数据集和logistic回归最佳拟合直线的函数
def plotBestFit(wei):
import matplotlib.pyplot as plt
weights=wei.getA() #将矩阵转换为数组
dataMat,labelMat=loadDataSet()
dataArr=array(dataMat)
n=shape(dataArr)[0]
xcord1=[]
ycord1=[]
xcord2=[]
ycord2=[]
for i in range(n):
if int(labelMat[i])==1:
xcord1.append(dataArr[i,1])
ycord1.append(dataArr[i,2])
else:
xcord2.append(dataArr[i,1])
ycord2.append(dataArr[i,2])
fig=plt.figure()
ax=fig.add_subplot(111)
ax.scatter(xcord1,ycord1,s=30,color='red',marker='s')
ax.scatter(xcord2, ycord2, s=30, color='green')
x=arange(-3,3,0.1)
y=(-weights[0]-weights[1]*x)/weights[2]
ax.plot(x,y)
plt.xlabel('x1')
plt.ylabel('x2')
plt.show()
运行程序：
if __name__=="__main__":
dataArr,labelMat=loadDataSet()
weights=gradAscent(dataArr,labelMat)
plotBestFit(weights)
运行结果：
可以看出虽然分类结果很好，但是尽管数据集很小，这个方法却需要大量的计算（300次乘法）。
5.2.4 训练算法：随机梯度上升
梯度上升法是每次更新回归系数时都需要遍历整个数据集，该方法在处理100个左右的数据及时尚可，但如果有数十亿样本和成千上万的特征，那么该方法的计算复杂度就太高了。一种改进方法是一次仅用一个样本点来更新回归系数，该方法称为随机梯度上升算法。由于可以在新样本到来时对分类器进行增量式更新，因而随机梯度上升法是一个在线学习算法。与“在线学习”相对应，一次处理所有数据被称作是“批处理”。
随机梯度上升法可以写成以下的伪代码：
所有回归系数初始化为1
对数据集中每个样本
。。计算该样本的梯度
。。使用alphaxgradient更新回归系数值
返回回归系数值
以下是随机梯度上升法的实现代码：
#随机梯度上升算法
def stocGradAscent0(dataMatrix,classLabels):
m,n=shape(dataMatrix)
alpha=0.01
weights=ones(n)
for i in range(m):
h=sigmoid(sum(dataMatrix[i]*weights))
error=classLabels[i]-h
weights=weights+alpha*error*dataMatrix[i]
return weights
可以看到，随机梯度上升法与梯度上升法在代码上很相似，但也有一些区别：第一，后者的变量h和误差error都是向量，而前者则全是数值；第二，前者没有矩阵的转换过程，所有变量的数据类型都是numpy数组。
为了验证该方法的结果，将上述代码添加到logRegres.py中，运行程序：
if __name__=="__main__":
dataArr,labelMat=loadDataSet()
weights=stocGradAscent0(array(dataArr),labelMat)
plotBestFit(weights)
运行结果：
可以看到随机梯度上升法拟合出来的结果不够完美，但是直接比较两个结果是不公平的，梯度上升法的结果是在整个数据集上迭代500次才得到的。一个判断优化算法优劣的可靠方法是看它是否收敛，也就是参数是否达到了稳定值，是否还会不断地变化？对此，在程序中对随机梯度上升法做一些修改，使其在整个数据集上运行150次。
#改进的随机梯度上升算法
def stocGradAscent1(dataMatrix,classLabels,numIter=150):
import random
m,n=shape(dataMatrix)
weights=ones(n)
for j in range(numIter):
dataIndex=range(m)
for i in range(m):
alpha=4/(1+j+i)+0.01
randIndex=int(random.uniform(0,len(dataIndex)))
h=sigmoid(sum(dataMatrix[randIndex]*weights))
error=classLabels[randIndex]-h
weights=weights+alpha*error*dataMatrix[randIndex]
return weights
在改进的算法中，‘alpha=4/(1+j+i)+0.01’使得alpha在每次迭代的时候都会调整，这回缓解回归系数波动或者高频波动。另外，虽然alpha会随着迭代次数不断减小，但永远不会减小到0，因为式子中还存在一个常数项。必须这样做的原因是为了保证在多次迭代之后新数据仍然具有一定的影响。如果要处理的问题是动态变化的，那么可以适当加大上述常数项，来确保新的值获更大的回归系数。另一点值得注意的是，在降低alpha的函数中，alpha每次减少1/(j+i),其中j是迭代次数，i是样本点的下标。这样当j<<max(i)时，alpha就不是严格下降的。避免参数的严格下降也常见于模拟退火算法等其它优化算法中。
在改进的程序中，通过随机选取样本来更新回归系数。这种方法将减少周期性的波动。这种方法每次随机从列表中选出一个值，然后从列表中删掉该值（再进行下一次迭代）。
运行结果：
可以看到随机梯度上升法迭代150次达到了和梯度上升算法相似的效果。
5.3 示例：从疝气病症预测病马的死亡率
数据集来源：https://archive.ics.uci.edu/ml/datasets/Horse+Colic：包含368个样本和28个特征，该数据集中包含了医院检测马疝病的一些指标，有的指标比较主观，有的指标难以测量，例如马的疼痛级别。
示例：使用logistic回归估计马疝病的死亡率
（1）收集数据：给定数据文件。
（2）准备数据：用python解析文本文件并填充缺失值。
（3）分析数据：可视化并观察数据。
（4）训练算法：使用优化算法，找到最佳的系数。
（5）测试算法：为了量化回归的效果，需要观察错误率。根据错误率决定是否回退到训练阶段，通过改变迭代的次数和步长等参数来得到更好的回归系数。
（6）使用算法：实现一个简单的命令行程序来收集马的症状并输出预测结果。
除了部分指标主观和难以测量外，该数据还存在一个问题，数据集中有30%的值是缺失的。下面将首先介绍如何处理数据集中的数据缺失问题，然后再利用logistic回归和随机梯度上升算法来预测病马的生死。
5.3.1 准备数据：处理数据中的缺失值
数据中的缺失值是个非常棘手的问题，有很多文献都致力于解决这个问题。那么，数据缺失究竟带来了什么问题？假设有100个样本和20个特征，这些数据都是机器学习收集回来的。若机器上的某个传感器损坏导致一个特征无效时该怎么办？此时是否要扔掉整个数据？这种情况下，另外19个特征怎么办？它们是否还可用？答案是肯定的。因为有时候数据相当昂贵，扔掉和重新获取都是不可取的，所以必须采用一些方法来解决这个问题。
下面是一些可选的做法：
*使用可用特征的均值来填补缺失值；
*使用特殊值来填补缺失值，如-1；
*忽略有缺失值的样本；
*使用相似样本的均值填补缺失值；
*使用另外的机器学习算法预测缺失值。
首先，需要对数据集进行预处理，使其可以顺利地使用分类算法。在预处理阶段需要做两件事：第一，所有缺失的值必须用一个实数值来替换，因为使用的numpy数据类型不允许包含缺失值。这里选择实数0来替换所有缺失值，恰好能适用于logistic回归。这样做不会在更新时不会影响系数的值。回归系数的更新公式如下：
如果dataMatrix的某特征对应值为0，那么该特征的系数将不做更新，即：
另外，由于sigmoid（0）=0.5，即它对结果的预测不具有任何倾向性，因此上述做法也不会对误差项造成任何影响。基于上述原因，缺失值用0代替可以保留现有数据，也不需要对优化算法进行修改。此外，数据集中的特征取值一般不为0，因此在某种意义上说它也满足“特殊值”这个要求。
预处理中做的第二件事是，如果在测试数据集中发现了一条数据的类别标签已经缺失，那么简单做法是将该条数据丢弃。这是因为类别标签与特征不同，很难确定采用某个合适的值来替换。采用logistic回归进行分类时这种做法是合理的，而如果采用类似kNN的方法就可能不太可行。
原始的数据集经过预处理之后保存成两个文件：horseColicTest.txt和horseColicTraining.txt。
5.3.2 测试算法：用logistic回归进行分类
使用logistic回归方法进行分类要做的是把测试集上每个特征向量乘以最优化方法得来的回归系数，再将该乘积结果求和，最后输入到sigmoid函数中即可。如果对应的sigmoid值大于0.5就预测类别标签为1，否则为0.
将下列代码添加到logRegres.py中：
def classifyVector(inX,weights):
prob=sigmoid(sum(inX*weights))
if prob>0.5:
return 1.0
else:
return 0.0
def colicTest():
frTrain=open('./data/horseColicTraining.txt')
frTest=open('./data/horseColicTest.txt')
trainingSet=[]
trainingLabels=[]
for line in frTrain.readlines():
currLine=line.strip().split(' ')
lineArr=[]
for i in range(27):
lineArr.append(float(currLine[i]))
trainingSet.append(lineArr)
trainingLabels.append(float(currLine[27]))
trainWeights=stocGradAscent1(array(trainingSet),trainingLabels,500)
errorCount=0
numTestVec=0.0
for line in frTest.readlines():
numTestVec+=1.0
currLine=line.strip().split(' ')
lineArr=[]
for i in range(27):
lineArr.append(float(currLine[i]))
if int(classifyVector(array(lineArr),trainWeights))!=int(currLine[27]):
errorCount+=1
errorRate=(float(errorCount)/numTestVec)
print('the error rate of this test is: %f'%errorRate)
return errorRate
def multiTest():
numTests=10
errorSum=0.0
for k in range(numTests):
errorSum+=colicTest()
print('after %d iterations the average error rate is: %f'%(numTests,float(errorSum)/float(numTests)))
第一个函数classifyVector（），它以回归系数和特征向量作为输入来计算对应的sigmoid值。如果sigmoid值大于0.5函数返回1，否则返回0.
接下来的函数colicTest（）是用于打开测试集和训练集，并对数据进行格式处理的函数。该函数首先导入训练集，最后一列是类别标签。数据最初有三个类别标签，分别代表马的三种情况：‘仍存活’、‘已经死亡’和‘已经安乐死’，这里为了方便，将‘已经死亡’和‘已经安乐死’合并成‘未能存活’这个标签。数据导入之后，便可以使用stocGradAscent1（）来计算回归系数向量。这里可以自由设定迭代的次数，例如在训练集上使用500次迭代，实验结果表明这比默认迭代150次的效果更好。在系数计算完成之后，导入测试集并计算分类错误率。整体看来，colicTest（）具有完全独立的功能，多次运行得到的结果可能稍有不同，这是因为其中有随机的成分在里面。如果在stocGradAscent1（）函数中回归系数已经完全收敛，那么结果才将是确定的。
最后一个函数是multiTest（），其功能是调用函数colicTest（）10次并求结果的平均值。
运行代码：
if __name__=="__main__":
multiTest()
运行结果（可以通过调整colicTest（）中的迭代次数和stocGradAscent（）中的步长改变平均错误率）：
the error rate of this test is: 0.632353
the error rate of this test is: 0.632353
the error rate of this test is: 0.632353
the error rate of this test is: 0.632353
the error rate of this test is: 0.632353
the error rate of this test is: 0.632353
the error rate of this test is: 0.632353
the error rate of this test is: 0.632353
the error rate of this test is: 0.632353
the error rate of this test is: 0.632353
after 10 iterations the average error rate is: 0.632353
※由于每次的结果都一样，所以可能有些错误，但是还未找到具体错误。
5.4 本章小结
logistic回归的目的是寻找一个非线性函数sigmoid（）的最佳拟合参数，求解过程可以由最优化算法来完成。在最优化算法中，最常用的是梯度上升算法，而梯度上升算法又可以简化为随机梯度上升算法。
随机梯度上升算法与梯度上升算法的效果相当，但占用更少的计算资源。此外，随机梯度上升是一个在线算法，它可以在新数据到来时就完成参数更新，而不需要重新读取整个数据集来进行批处理运算。
机器学习的一个重要问题就是如何处理缺失数据。这个问题没有标准答案，取决于实际应用中的需求。现有一些解决方案，每种方案都各有优缺点。