目录
3、Logistic回归算法实例1--从疝气病预测病马的死亡率
1、Logistic回归
本篇首先阐述Logistic回归的定义，然后介绍一些最优化算法，其中包括基本的梯度上升法和一个改进的随机梯度上升法，这些最优化算法将用于分类器的训练。本章最后会给出一个Logistic回归的实例，预测一匹病马是否能被治愈。
我们想要的函数应该是，能接受所有的输入然后预测出类别。例如，在两个类的情况下，上述函数输出0或1。或许你之前接触过具有这种性质的函数，该函数称为海维塞德阶跃函数( Heaviside step function )，或者直接称为单位阶跃函数。然而，海维塞德阶跃函数的问题在于:该函数在跳跃点上从0瞬间跳跃到1，这个瞬间跳跃过程有时很难处理。幸好，另一个函数也有类似的性质R，且数学上更易处理，这就是Sigmoid函数。Sigmoid函数具体的计算公式如下:
图下给出了Sigmoid函数在不同坐标尺度下的两条曲线图。当x为0时，Sigmoid函数值为0.5。随着x的增大，对应的Sigmoid值将逼近于1;而随着x的减小，Sigmoid值将逼近于0。如果横坐标刻度足够大，Sigmoid函数看起来很像一个阶跃函数。
因此，为了实现Logistic回归分类器，我们可以在每个特征上都乘以一个回归系数，然后把所有的结果值相加，将这个总和代人Sigmoid函数中，进而得到-一个范围在0~1之间的数值。任何大于0.5的数据被分入1类，小于0.5即被归入0类。所以，Logistic回归也可以被看成是一种概率估计。
确定了分类器的函数形式之后，现在的问题变成了:最佳回归系数是多少?如何确定它们的大小?
Sigmoid函数的输入记为z，由下面公式得出:
如果采用向量的写法，上述公式可以写成，它表示将这两个数值向量对应元素相乘然后全部加起来即得到z值。其中的向量x是分类器的输入数据，向量w也就是我们要找到的最佳参数(系数)从而使得分类器尽可能地精确。为了寻找该最佳参数,需要用到最优化理论的一些知识。
下面首先介绍梯度上升的最优化方法，我们将学习到如何使用该方法求得数据集的最佳参数。接下来，展示如何绘制梯度上升法产生的决策边界图，该图能将梯度上升法的分类效果可视化地呈现出来。最后我们将学习随机梯度上升算法，以及如何对其进行修改以获得更好的结果。
你最经常听到的应该是梯度下降算法，它与这里的梯度上升算法是一样的，只是公式中的加法需要变成减法。
梯度上升对应的公式可以写成：
梯度下降对应的公式可以写成：
他们的区别是梯度上升算法用来求函数的最大值，而梯度下降算法用来求函数的最小值。
梯度上升方法的原理可以参考梯度下降：机械学习的训练模型（线性回归、梯度下降、正规方程）_吃花生的图图的博客-CSDN博客
2、Logistic回归代码
基于上面的内容，我们来看一个Logistic回归分类器的应用例子.
图5-3中有100个样本点，每个点包含两个数值型特征:X1和X2。在此数据集上，我们将通过使用梯度上升法找到最佳回归系数，也就是拟合出Logistic回归模型的最佳参数。
1、预处理数据
数据在这：数据在这！ 提取码：9hip
#数据加载
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,labelMat2、sigmoid函数
#Sigmoid函数
def sigmoid(inX):
return 1.0/(1+np.exp(-inX))3、回归梯度上升函数
变量alpha是向目标移动的步长，maxCycles是迭代次数。在for循环迭代完成后，将返回训练好的回归系数。需要强调的是，运算是矩阵运算。变量h不是一个数而是一个列向量，列向量的元素个数等于样本个数，这里是100。对应地，运算dataMatrix * weights代表的不止一次乘积计算，事实上该运算包含了300次的乘积。
def gradAscent(dataMathIn,classLabel):
dataMatrix = np.mat(dataMathIn) #将列表转化为numpy的矩阵
labelMat = np.mat(classLabel).transpose() #将列表转化为矩阵，并求转置（transpose)
m,n = np.shape(dataMatrix)
alpha = 0.001
maxCycles = 500 #迭代次数
weight = np.ones((n,1))
for k in range(maxCycles):
h = sigmoid(dataMatrix*weight)
error = labelMat - h #(y-h)误差值
weight = weight + alpha*dataMatrix.transpose()*error #w=w+alpha*(y-h)*x梯度上升的迭代公式
return weight4、画出决策边界
#画出决策边界
def plotBestFit(weight):
import matplotlib.pyplot as plt
dataMat,labelMat = loadDataSet()
dataArr = np.array(dataMat)
n = np.shape(dataArr)[0] #[0]求dataArr的行数，1，是求列数
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) #添加subplot
ax.scatter(xcord1,ycord1,s = 30,c='red',marker = 's') #绘制正样本
ax.scatter(xcord2,ycord2,s = 30, c = 'green') #绘制负样本
w = np.array(weight)
x = np.arange(-3.0,3.0,0.1) #X轴的取值范围
y = (-w[0] - w[1] * x) / w[2] #设sigmoid =0,因为sigmoid = 0是两种类别的分界线 sigmoid = w0x0+w1x1+w2x2 = 0
ax.plot(x,y)
plt.xlabel('X1')
plt.ylabel('X2')
plt.show()画出的图像为：
这个分类结果相当不错,从图上看只错分了两到四个点。但是,尽管例子简单且数据集很小，这个方法却需要大量的计算（300次乘法)。因此下一节将对该算法稍作改进，从而使它可以用在真实数据集上。
5、随机梯度上升
梯度上升算法在每次更新回归系数时都需要遍历整个数据集，该方法在处理100个左右的数据集时尚可，但如果有数十亿样本和成千上万的特征，那么该方法的计算复杂度就太高了。一种改进方法是一次仅用一个样本点来更新回归系数，该方法称为随机梯度上升算法。由于可以在新样本到来时对分类器进行增量式更新，因而随机梯度上升算法是一个在线学习算法。与“在线学习”相对应，一次处理所有数据被称作是“批处理”。
具体的原理可以查看：随机梯度下降、Map-reduce_吃花生的图图的博客-CSDN博客
#随机梯度上升
def stoGradAscent0(dataMatrix,classLabels):
dataMatrix=np.array(dataMatrix)
m,n = np.shape(dataMatrix)
alpha = 0.01
weight = np.ones(n)
for i in range(m):
h = sigmoid(sum(dataMatrix[i]*weight))
error = classLabels[i] - h
weight = weight +alpha*error*dataMatrix[i]
return weight执行完毕后将得到图下所示的最佳拟合直线图，该图与之前的回归梯度上升函数绘制出来的图有一些相似之处。可以看到，拟合出来的直线效果还不错，但并不像图5-4那样完美。这里的分类器错分了三分之一的样本。
直接比较两段代码结果是不公平的，后者的结果是在整个数据集上迭代了500次才得到的。一个判断优化算法优劣的可靠方法是看它是否收敛，也就是说参数是否达到了稳定值，是否还会不断地变化?对此，我们在程序清单5-3中随机梯度上升算法上做了些修改，使其在整个数据集上运行200次。最终绘制的三个回归系数的变化情况如图下所示。
值得注意的是，在大的波动停止后，还有一些小的周期性波动。不难理解，产生这种现象的原因是存在一些不能正确分类的样本点(数据集并非线性可分)，在每次迭代时会引发系数的剧烈改变。我们期望算法能避免来回波动，从而收敛到某个值。另外，收敛速度也需要加快。
6、改进的随机梯度上升
第一处改进，一方面，alpha在每次迭代的时候都会调整，这会缓解图上的数据波第一处改进在①处。一方面，alpha在每次迭代的时候都会调整，这会缓解图上的数据波动或者高频波动。另外,虽然alpha会随着迭代次数不断减小，但永远不会减小到0,这是因为还存在一个常数项。必须这样做的原因是为了保证在多次迭代之后新数据仍然具有一定的影响。如果要处理的问题是动态变化的，那么可以适当加大上述常数项，来确保新的值获得更大的回归系数。另一点值得注意的是，在降低alpha的函数中，alpha每次减少1/(j+i),其中j是迭代次数，i是样本点的下标。这样当j<<max(i)时，alpha就不是严格下降的。避免参数的严格下降也常见于模拟退火算法等其他优化算法中。
第二个改进的地方，这里通过随机选取样本来更新回归系数。这种方法将减少周期性的波动。具体实现方法与第3章类似，这种方法每次随机从列表中选出一个值，然后从列表中删掉该值（再进行下一次迭代)。
#随机梯度算法的改进
def stoGradAscent1(dataMatrix,classLabels,numIter = 150):
dataMatrix=np.array(dataMatrix)
m,n = np.shape(dataMatrix)
weight = 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 #降低alpha的大小，每次减小1/(j+i)。
randIndex = int(random.uniform(0,len(dataIndex))) #随机选取样本,可以实现随机选取更新
h = sigmoid(sum(dataMatrix[randIndex]*weight))
error = classLabels[randIndex] - h
weight = weight + alpha*error*dataMatrix[randIndex]
del(dataIndex[randIndex])
return weight程序运行之后应该能看到类似图下结果：
7、完整代码
import numpy as np
import array
import random
#数据加载
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
#Sigmoid函数
def sigmoid(inX):
return 1.0/(1+np.exp(-inX))
def gradAscent(dataMathIn,classLabel):
dataMatrix = np.mat(dataMathIn) #将列表转化为numpy的矩阵
labelMat = np.mat(classLabel).transpose() #将列表转化为矩阵，并求转置（transpose)
m,n = np.shape(dataMatrix)
alpha = 0.001
maxCycles = 500 #迭代次数
weight = np.ones((n,1))
for k in range(maxCycles):
h = sigmoid(dataMatrix*weight)
error = labelMat - h #(y-h)误差值
weight = weight + alpha*dataMatrix.transpose()*error #w=w+alpha*(y-h)*x梯度上升的迭代公式
return weight
#画出决策边界
def plotBestFit(weight):
import matplotlib.pyplot as plt
dataMat,labelMat = loadDataSet()
dataArr = np.array(dataMat)
n = np.shape(dataArr)[0] #[0]求dataArr的行数，1，是求列数
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) #添加subplot
ax.scatter(xcord1,ycord1,s = 30,c='red',marker = 's') #绘制正样本
ax.scatter(xcord2,ycord2,s = 30, c = 'green') #绘制负样本
w = np.array(weight)
x = np.arange(-3.0,3.0,0.1) #X轴的取值范围
y = (-w[0] - w[1] * x) / w[2] #设sigmoid =0,因为sigmoid = 0是两种类别的分界线 sigmoid = w0x0+w1x1+w2x2 = 0
ax.plot(x,y)
plt.xlabel('X1')
plt.ylabel('X2')
plt.show()
#随机梯度上升
def stoGradAscent0(dataMatrix,classLabels):
dataMatrix=np.array(dataMatrix)
m,n = np.shape(dataMatrix)
alpha = 0.01
weight = np.ones(n)
for i in range(m):
h = sigmoid(sum(dataMatrix[i]*weight))
error = classLabels[i] - h
weight = weight +alpha*error*dataMatrix[i]
return weight
#随机梯度算法的改进
def stoGradAscent1(dataMatrix,classLabels,numIter = 150):
dataMatrix=np.array(dataMatrix)
m,n = np.shape(dataMatrix)
weight = 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 #降低alpha的大小，每次减小1/(j+i)。
randIndex = int(random.uniform(0,len(dataIndex))) #随机选取样本,可以实现随机选取更新
h = sigmoid(sum(dataMatrix[randIndex]*weight))
error = classLabels[randIndex] - h
weight = weight + alpha*error*dataMatrix[randIndex]
del(dataIndex[randIndex])
return weight
dataMat,labelMat = loadDataSet()
# weight=gradAscent(dataMat,labelMat)
weight = stoGradAscent1(dataMat,labelMat)
plotBestFit(weight)3、Logistic回归算法实例1--从疝气病预测病马的死亡率
数据：数据在这！ 提取码：8z35
前面几节介绍了优化算法，但目前为止还没有在分类上做任何实际尝试。使用Logistic回归方法进行分类并不需要做很多工作,所需做的只是把测试集上每个特征向量乘以最优化方法得来的回归系数，再将该乘积结果求和，最后输入到Sigmoid函数中即可。如果对应的Sigmoid值大于0.5就预测类别标签为1，否则为0。
#分类器
def classifyVector(inX,weight):
inX = inX.astype(np.float64)
weight = np.array(weight)
t = sum(inX*weight)
prob = sigmoid(t)
if prob>0.5:
return 1.0
else:
return 0.0
def 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]))
trainingWeight =stoGradAscent1(np.array(trainingSet),trainingLabels,500)
errorCount = 0
numTestVec = 0.0
for line in frTest.readlines():
numTestVec +=1.0
currLine = line.strip().split('\t')
lineArr = []
for i in range(21):
lineArr.append(currLine[i])
if (int(classifyVector(np.array(lineArr),trainingWeight))!=int(currLine[21])):
errorCount += 1
errorRate = float(errorCount)/numTestVec
print('The error rate of this test is:%f'%errorRate)
return errorRate
#调用10次colicTest()，求算法错误率的平均值
def mutiTest():
numTests = 10
errorsum = 0.0
for k in range(numTests):
errorsum +=colicTest()
print("after %d iteration the average error is:%f"%numTests,(errorsum/float(numTests))*100)
从上面的结果可以看到，10次迭代之后的平均错误率为35%。事实上，这个结果并不差，因为有30%的数据缺失。当然，如果调整colicTest()中的迭代次数和stochGradAscent1()中的步长，平均错误率可以降到20%左右。
4、小结
Logistic回归的目的是寻找一个非线性函数Sigmoid的最佳拟合参数，求解过程可以由最优化算法来完成。在最优化算法中，最常用的就是梯度上升算法，而梯度上升算法又可以简化为随机梯度上升算法。
随机梯度上升算法与梯度上升算法的效果相当，但占用更少的计算资源。此外，随机梯度上升是一个在线算法，它可以在新数据到来时就完成参数更新，而不需要重新读取整个数据集来进行批处理运算。
4956
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
