python机器学习实战（四）
版权声明：本文为博主原创文章，转载请指明转载地址
http://www.cnblogs.com/fydeblog/p/7364317.html
前言
这篇notebook是关于机器学习中logistic回归，内容包括基于logistic回归和sigmoid分类，基于最优化方法的最佳系数确定，从疝气病症预测病马的死亡率。操作系统：ubuntu14.04 运行环境：anaconda-python2.7-jupyter notebook 参考书籍：机器学习实战和源码 notebook writer —-方阳
注意事项：在这里说一句，默认环境python2.7的notebook，用python3.6的会出问题，还有我的目录可能跟你们的不一样，你们自己跑的时候记得改目录，我会把notebook和代码以及数据集放到结尾的百度云盘，方便你们下载！
1. 基于logistic回归和sigmoid函数的分类
首先说说sigmoid函数吧
它的表达式是 g(z) = 1/(1+exp(-x)) ，为直观看出，我们画画这个函数的曲线
import matplotlib.pyplot as plt
import numpy as np x = np.linspace(-5,5,200)
y = 1./(1+np.exp(-x)) plt.figure()
plt.plot(x,y)
plt.xlabel('x')
plt.ylabel('sigmiod(x)')
plt.show()
上面就是sigmiod函数的图形，那么我们怎么用sigmiod函数进行逻辑回归判决呢？
首先观察函数图形，sigmiod函数的y轴被限制在区间（0，1）上，这有利于我们判决，将线性的无穷范围压缩到这个小范围，当x=0的时候，sigmiod(0) = 0.5, 于是我们就将0.5当作界限，特征值乘以一个回归系数，然后结果相加，代入到这个sigmiod函数当中，将函数值大于0.5分为1类，小于0.5的分为0类，至此，logistic分类完成。
2. 基于最优化方法的最佳回归系数确定
这里说一下，sigmiod函数是为了帮助我们来判断分类类别，然后与真实类别相比较，算出误差，然后用梯度上升最小化误差，得到最佳系数。
sigmiod函数的输入记为z，公式：z = w0x0+w1x1+w2x2+…+wnxn (这里0，1，2，…,n都代表下标系数)，简单写就是 z = wTx （T代表转置）
2.1 梯度上升法
梯度上升法的思想：要找到某函数的最大值 ,最好的方法是沿着该函数的梯度方向探寻。
这是函数f(x,y)的梯度表达式，当要沿x方向移动时，就是对x求偏导；当要沿y方向移动时，就是对y求偏导。其中，函数f(x,y)必须要在待计算的点上可微。
梯度上升算法的迭代公式：w := w+a dw（f(w)） (dw是关于系数w的梯度，a是学习率)
梯度上升法的具体上升过程，如图所示
该公式将一直被迭代执行,直至达到某个停止条件为止,比如迭代次数达到某个指定值或算法达到某个可以允许的误差范围。
梯度上升算法与梯度下降算法的区别：就是梯度上升的学习率前面的加号变为减号就是梯度下降算法
2.2 训练算法 :使用梯度上升找到最佳参数
梯度上升法的伪代码如下:
每个回归系数初始化为1
重复R次：
计算整个数据集的梯度
使用alpha X gradient更新回归系数
返回回归系数
下面进入梯度上升法的具体实现
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) #convert to NumPy matrix
labelMat = mat(classLabels).transpose() #convert to NumPy matrix
m,n = shape(dataMatrix)
alpha = 0.001
maxCycles = 500
weights = ones((n,1))
for k in range(maxCycles): #heavy on matrix operations
h = sigmoid(dataMatrix*weights) #matrix mult
error = (labelMat - h) #vector subtraction
weights = weights + alpha * dataMatrix.transpose()* error #matrix mult
return weights
第一个函数loadDataSet函数就是导入数据集，并进行封装，将特征值封装成三列的多维列表dataMat，label放在labelMat里面
第二个函数不用多说，就是sigmiod函数
第三个函数就是梯度上升的具体算法，dataMatIn就是上面的多维列表dataMat，classlabels就是labelMat，函数首先将输入的数据集和标签全部转化成的numpy矩阵，是为了能够进行矩阵运算和向量运算，maxCycles代表迭代的最大次数，至于参数的迭代还可以看看以下图片，我选的是吴恩达的ppt上，这写的是梯度下降的迭代过程（梯度上升就是反过来，原理一样），可见损失函数的梯度会有逐渐趋于0，这样迭代公式的梯度那项也会趋于0，参数不会有太大的浮动，趋于稳定，误差最小。
测试一下吧
cd 桌面/machinelearninginaction/Ch05
/home/fangyang/桌面/machinelearninginaction/Ch05
import logRegres
dataMat , labelMat = logRegres.loadDataSet()
logRegres.gradAscent(dataMat,labelMat)
2.3 分析数据：画出决策边界
我们得到了最佳系数，便于理解，我们要画出分隔线，便于我们观察
代码如下：
def plotBestFit(weights):
import matplotlib.pyplot as plt
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, c='red', marker='s')
ax.scatter(xcord2, ycord2, s=30, c='green')
x = arange(-3.0, 3.0, 0.1)
y = (-weights[0]-weights[1]*x)/weights[2]
ax.plot(x, y)
plt.xlabel('X1'); plt.ylabel('X2');
plt.show()
这个函数先是调用loadDataSet函数将数据集和标签赋给dataMat,labelMat，然后对不同类别进行不同的分组，类别1的数据放在xcord1和ycord1，类别2的数据放在xcord2和ycord2，然后分别显示，最后画出输入的权重对应的分隔线，y的求解你可能有疑问，这里说一下，具体表达式是wTx=0,wT是输入权重，x=[x0,x1,x2],其中x0为了方便表示所建立的，值为1，x1就是上述函数的x，x2就是y，这下你就知道为什么是那个表达式了吧。
from numpy import *
weights = logRegres.gradAscent(dataMat,labelMat)
logRegres.plotBestFit(weights.getA()) # the function of getA is used to transform matrix into array
可以看出图中只错分了两到四个点，效果不错。
2.4 训练算法：随机梯度上升
梯度上升算法在每次更新回归系数时都需要遍历整个数据集，当特征的数目非常多的时候，计算量会非常巨大。
一种改进方法是一次仅用一个样本点来更新回归系数, 该方法称为随机梯度上升算法。
由于可以在新样本到来时对分类器进行增量式更新,因而随机梯度上升算法是一个在线学习算法。与在线学习相对应 ,一次处理所有数据被称作是 “批处理” 。
随机梯度上升算法的伪代码如下
所有回归系数初始化为 1
对数据集中每个样本
计算该样本的梯度
使用 alpha x gradient更新回归系数值
返回回归系数值
代码如下：
def stocGradAscent0(dataMatrix, classLabels):
m,n = shape(dataMatrix)
alpha = 0.01
weights = ones(n) #initialize to all ones
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数组。
这是因为梯度上升算法是遍历所有数据，形成的是所有数据的向量，而随机梯度上升算法每次只用一个样本，所以是单一数值。其他类似
来测试一下效果吧
weights = logRegres.stocGradAscent0(array(dataMat),labelMat)
logRegres.plotBestFit(weights)
可见拟合的直线不完美，错分了三分之一的样本，直接比较结果是不公平的，梯度上升算法是在整个数据集中迭代了500次，而随机梯度算法只在整个数据集中迭代了1次，计算量相差很多倍。所以这里还需对随机梯度上升算法进行优化，代码如下
def stocGradAscent1(dataMatrix, classLabels, numIter=150):
m,n = shape(dataMatrix)
weights = ones(n) #initialize to all ones
for j in range(numIter):
dataIndex = range(m)
for i in range(m):
alpha = 4/(1.0+j+i)+0.0001 #apha decreases with iteration, does not
randIndex = int(random.uniform(0,len(dataIndex)))
h = sigmoid(sum(dataMatrix[randIndex]*weights))
error = classLabels[randIndex] - h
weights = weights + alpha * error * dataMatrix[randIndex]
del(dataIndex[randIndex])
return weights
这里主要改了两点：
第一点增加了alpha动态减少的机制，这样做的原因是为了保证在多次迭代之后新数据仍然具有一定的影响。
第二点是通过随机选取样本来更新回归系数，这种方法减少周期性的波动。每次随机从列表中选出一个值,然后从列表中删掉该值，重新迭代
需要注意的是：
如果要处理的问题是动态变化的,那么可以适当加大上述常数项,来确保新的值获得更大的回归系数。
再次运行出来看看
weights = logRegres.stocGradAscent1(array(dataMat),labelMat)
logRegres.plotBestFit(weights)
可见分类效果与梯度上升算法差不多，比较一下计算量，梯度上升算法迭代了500次的整个数据集，而随机梯度上升算法迭代了150次就达到类似的效果。
3. 示例:从疝气病症预测病马的死亡率
这个例子是通过马疝病的一些指标，使用logistic回归和随机梯度上升算法来预测病马的生死。
3.1 准备数据：处理数据中的缺失值
马疝病的数据集中有30%的值是缺失的，我们怎样来解决这个问题呢？
首先我们要知道，有时候数据是非常昂贵的，扔掉缺失数据和重新获取新的数据都是不可取的，所以我们采用一些方法来解决这个问题，方法如下：
下面给出了一些可选的做法:
□ 使用可用特征的均值来填补缺失值;
□ 使用特殊值来补缺失值,如 -1;
□ 忽略有缺失值的样本;
□ 使用相似样本的均值添补缺失值;
□ 使用另外的机器学习算法预测缺失值。
预处理阶段的两件事：
第一件事，所有的缺失值必须用一个实数值来替换,因为我们使用的numpy数据类型不允许包含缺失值。
第二件事，如果在测试数据集中发现了一条数据的类别标签已经缺失,那么我们的简单做法是将该条数据丢弃。这是因为类别标签与特征不同,很难确定采用某个合适的值来替换。
这个例子选实数0来替换所有缺失值，不影响特征系数，如果等于0，对应的参数也被置0，不会更新，还有就是sigmiod(0)=0.5，即对结果的预测不具有任何倾向性，所有用0代替缺失值
3.2 测试算法 :用Logistic回归进行分类
使用Logistic回归方法进行分类并不需要做很多工作,所需做的只是把测试集上每个特征向量乘以最优化方法得来的回归系数,再将该乘积结果求和,最后输人到sigmiod函数中即可。如果对应的sigmiod值大于0.5就预测类别标签为1 , 否则为0 。
例子的代码如下
def classifyVector(inX, weights):
prob = sigmoid(sum(inX*weights))
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]))
trainWeights = stocGradAscent1(array(trainingSet), trainingLabels, 1000)
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(float(currLine[i]))
if int(classifyVector(array(lineArr), trainWeights))!= int(currLine[21]):
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, errorSum/float(numTests))
第一个函数classifyVector，它以回归系数和特征向量作为输入来计算对应的sigmiod值，如果值大于0.5返回1，否则返回0
第二个函数colicTest，先是打开训练集和测试集，然后使用改善后的随机梯度上升算法对训练集进行训练，得到训练参数，然后使用这个训练参数带入到测试集当中，比较这个测试集使用该参数得到的类别和实际类别，算出错误的个数，最终程序是返回错误率
第三个函数就是多次使用第二个函数，最后打印出平均错误率
logRegres.multiTest()
平均错误率34.7761%，这个错误率还好，毕竟我们有30%的数据缺失嘛，作者书上调整迭代次数和步长，可以改善错误率，你们可以试一试。
小结
1. logistic回归的目的是寻找一个非线性函数sigmiod的最佳拟合参数,求解过程可以由最优化算法来完成。在最优化算法中,最常用的就是梯度上升算法, 而梯度上升算法又可以简化为随机梯度上升算法。
2. 随机梯度上升算法与梯度上升算法的效果相当, 但占用更少的计算资源。此 外 ,随机梯度上升是一个在线算法, 它可以在新数据到来时就完成参数更新, 而不需要重新读取整个数据集来进行批处理运算。
百度云链接：https://pan.baidu.com/s/1cSRMHwvic0ujKL-oWqQ1aA
