logistic回归
一、logistic回归概述
主要是进行二分类预测，也即是对于0~1之间的概率值，当概率大于0.5预测为1，小于0.5预测为0.显然，我们不能不提到一个函数，即sigmoid=1/(1+exp(-inX)),该函数的曲线类似于一个s型，在x=0处，函数值为0.5.
于是，为了实现logistic分类器，我们可以在每个特征上都乘以一个回归系数，然后所有的相乘结果进行累加，将这个总结作为输入，输入到sigmoid函数中，从而得到一个大小为0~1之间的值，当该值大于0.5归类为1，否则归类为0，这样就完成了二分类的任务。所以logistic回归可以看成是一种概率估计。
二，基于最优化方法的最佳回归系数确定
sigmoid函数的输入记为z，即z=w0x0+w1x1+w2x2+...+wnxn，如果用向量表示即为z=wTx，它表示将这两个数值向量对应元素相乘然后累加起来。其中，向量x是分类器的输入数据，w即为我们要拟合的最佳参数，从而使分类器预测更加准确。也就是说，logistic回归最重要的是要找到最佳的拟合参数。
1 梯度上升法
梯度上升法(等同于我们熟知的梯度下降法，前者是寻找最大值，后者寻找最小值)，它的基本思想是：要找到某函数的最大值，最好的方法就是沿着该函数的梯度方向搜寻。如果函数为f，梯度记为D，a为步长，那么梯度上升法的迭代公式为：w：w+a*Dwf(w)。该公式停止的条件是迭代次数达到某个指定值或者算法达到某个允许的误差范围。我们熟知的梯度下降法迭代公式为：w：w-a*Dwf(w)
2 使用梯度上升法寻找最佳参数
假设有100个样本点，每个样本有两个特征：x1和x2.此外为方便考虑，我们额外添加一个x0=1，将线性函数z=wTx+b转为z=wTx(此时向量w和x的维度均价1).那么梯度上升法的伪代码如下：
初始化每个回归系数为1
重复R次：
计算整个数据集梯度
使用alpha*gradient更新回归系数的向量
返回回归系数
代码如下：
#预处理数据 def loadDataSet(): dataMat=[];labelMat=[] fr=open('testSet.txt') #遍历文本的每一行 for line in fr.readlines(): #对当前行去除首尾空格，并按空格进行分离 lineArr=line.strip().split() #将每一行的两个特征x1，x2，加上x0=1,组成列表并添加到数据集列表中 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+exp(-inx)) #梯度上升法更新最优拟合参数 #@dataMatIn：数据集 #@classLabels：数据标签 def gradAscent(dataMatIn,classLabels): #将数据集列表转为Numpy矩阵 dataMatrix=mat(dataMatIn) 将数据集标签列表转为Numpy矩阵，并转置 labelMat=mat(classLabels).transpose() #获取数据集矩阵的行数和列数 m,n=shape(dataMat) #学习步长 alpha=0.001 #最大迭代次数 maxCycles=500 #初始化权值参数向量每个维度均为1.0 weights=one((n,1)) #循环迭代次数 for k in range(maxCycles): #求当前的sigmoid函数预测概率 h=sigmoid(dataMatrix*weights) #*********************************************** #此处计算真实类别和预测类别的差值 #对logistic回归函数的对数释然函数的参数项求偏导 error=(labelMat-h) #更新权值参数（这里的梯度是由最大似然损失函数的梯度求得） weights=weights+alpha*dataMatrix.transpose()*error #*********************************************** return weights
我们知道对回归系数进行更新的公式为：w：w+alpha*gradient,其中gradient是对参数w求偏导数。则我们可以通过求导验证logistic回归函数对参数w的梯度为(yi-sigmoid(wTx))*x,证明过程如下所示：
我们还可以通过matpotlib画出决策的边界，集体代码为：
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, 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()
3 随机梯度上升法
我们知道梯度上升法每次更新回归系数都需要遍历整个数据集，当样本数量较小时，该方法尚可，但是当样本数据集非常大且特征非常多时，那么随机梯度下降法的计算复杂度就会特别高。一种改进的方法是一次仅用一个样本点来更新回归系数，即岁集梯度上升法。由于可以在新样本到来时对分类器进行增量式更新，因此随机梯度上升法是一个在线学习算法。随机梯度上升法可以写成如下伪代码：
所有回归系数初始化为1
对数据集每个样本
计算该样本的梯度
使用alpha*gradient更新回顾系数值
返回回归系数值
#梯度上升算法 def stocGradAscent(dataMatrix,classLabels): #为便于计算，转为Numpy数组 dataMat=array(dataMatrix) 获取数据集的行数和列数 m,n=shape(dataMatrix) #初始化权值向量各个参数为1.0 weights=ones(n) #设置步长为0.01 alpha=0.01 #循环m次，每次选取数据集一个样本更新参数 for i in range(m): #计算当前样本的sigmoid函数值 h=sigmoid(dataMatrix[i]+weights) #计算当前样本的残差(代替梯度) error=(classLabels[i]-h) #更新权值参数 weights=weights+alpha*error*dataMatrix[i] return weights
评判一个优化算法的优劣的可靠方法是看其是否收敛，也就是说参数的值是否达到稳定值。此外，当参数值接近稳定时，仍然可能会出现一些小的周期性的波动。这种情况发生的原因是样本集中存在一些不能正确分类的样本点(数据集并非线性可分)，所以这些点在每次迭代时会引发系数的剧烈改变，造成周期性的波动。显然我们希望算法能够避免来回波动，从而收敛到某个值，并且收敛速度也要足够快。
为此，需要对上述随机梯度上升法代码进行适当修改，代码如下：
#@dataMatrix：数据集列表 #@classLabels：标签列表 #@numIter：迭代次数，默认150 def stocGradAscent1(dataMatrix,classLabels,numIter=150): #将数据集列表转为Numpy数组 dataMat=array(dataMatrix) #获取数据集的行数和列数 m,n=shape(dataMat) #初始化权值参数向量每个维度均为1 weights=ones(n) #循环每次迭代次数 for j in range(numIter): #获取数据集行下标列表 dataIndex=range(m) #遍历行列表 for i in range(m): #每次更新参数时设置动态的步长，且为保证多次迭代后对新数据仍然具有一定影响 #添加了固定步长0.1 alpha=4/(1.0+j+i)+0.1 #随机获取样本 randomIndex=int(random.nuiform(0,len(dataIndex))) #计算当前sigmoid函数值 h=sigmoid(dataMat[randomIndex]*weights) #计算权值更新 #*********************************************** error=classLabels-h weights=weights+alpha*error*dataMat[randomIndex] #*********************************************** #选取该样本后，将该样本下标删除，确保每次迭代时只使用一次 del(dataIndex[randomIndex]) return weights
上述代码中有两处改进的地方：
1 alpha在每次迭代更新是都会调整，这会缓解数据波动或者高频运动。此外，alpha还有一个常数项，目的是为了保证在多次迭代后仍然对新数据具有一定的影响，如果要处理的问题是动态变化的，可以适当加大该常数项，从而确保新的值获得更大的回归系数。
2 第二个改进的地方是选择随机的样本对参数进行更新，由于增加了随机性，这就防止参数值发生周期性的波动。
3 并且采用上述改进算法后，收敛速度更快
三，logistic回归实例：从疝气病症预测病马死亡率
现有数据集有100个样本和20个特征，但是数据存在一定的问题，即数据集有30%的缺失，因此，我们在对病马进行预测死亡率前，首先要解决数据的缺失问题。
我们可能会遇到数据缺失的情况，但有时候数据相当昂贵，扔掉和重新获取均不可取，这显然是会带来更多的成本负担，所以我们可以选取一些有效的方法来解决该类问题。比如：
1 使用可用特征的均值填补缺失值
2 使用特殊值来填补缺失的特征，如-1
3 忽略有缺失值的样本
4 使用相似样本的平均值填补缺失值
5 使用另外的机器学习算法预测缺失值
这里我们根据logstic回归的函数特征，选择实数0来替换所有缺失值，而这恰好能适用logistic回归。因此，它在参数更新时不会影响参数的值。即如果某特征对应值为0 ，那么由公式w：w+alpha*gradient，可知w不会发生改变。
此外，由于sigmoid(0)=0.5,表面该特征对结果的预测不具有任何倾向性，因此不会对误差造成影响。
当然，如果是发生有样本的类标签缺失的情况，此时我们最好的办法是将该样本舍弃，这是因为标签与特征不同，我们很难确定采用某个合适的值替换掉。
下面来看回归分类器代码：
#------------------------------实例：从疝气病预测病马的死亡率---------------------------- #1 准备数据：处理数据的缺失值 #这里将特征的缺失值补0，从而在更新时不影响系数的值 #2 分类决策函数 def clasifyVector(inX,weights): #计算logistic回归预测概率 prob=sigmoid(inX*weights) #大于0.5预测为1 if prob>0.5: return 1.0 #否则预测为0 else: return 0.0 #logistic回归预测算法 def colicTest(): #打开训练数据集 frTrain=open('horseColicTraining.txt') #打开测试数据集 frTest=open('horseColicTest.txt') #新建两个孔列表，用于保存训练数据集和标签 trainingSet=[];trainingLabels=[] #读取训练集文档的每一行 for line in frTrain.readlines(): #对当前行进行特征分割 currLine=line.strip().split() #新建列表存储每个样本的特征向量 lineArr=[] #遍历每个样本的特征 for i in range(21): 将该样本的特征存入lineArr列表 lineArr.append(float(currLine[i])) #将该样本标签存入标签列表 trainingLabels.append(currLine[21]) #将该样本的特征向量添加到数据集列表 trainingSet.append(lineArr) #调用随机梯度上升法更新logistic回归的权值参数 trainWeights=stocGradAscent1(trainingSet,trainingLabels,500) #统计测试数据集预测错误样本数量和样本总数 errorCount=0;numTestVec=0.0 #遍历测试数据集的每个样本 for line in frTest.readlines(): #样本总数加1 numTestVec+=1.0 #对当前行进行处理，分割出各个特征及样本标签 currLine=line.strip().split() #新建特征向量 lineArr=[] #将各个特征构成特征向量 for i in range(21): lineArr.append(float(currLine[i])) #利用分类预测函数对该样本进行预测，并与样本标签进行比较 if(clasifyVector(lineArr,trainWeights)!=currLine[21]): #如果预测错误，错误数加1 errorCount+=1 #计算测试集总的预测错误率 errorRate=(float(errorCount)/numTestVec) #打印错误率大小 print('the error rate of this test is: %f', %(errorRate)) #返回错误率 return errorRate #多次测试算法求取预测误差平均值 def multTest(): #设置测试次数为10次，并统计错误率总和 numTests=10;errorRateSum=0.0 #每一次测试算法并统计错误率 for k in range(numTests): errorRateSum+=colicTest() #打印出测试10次预测错误率平均值 print('after %d iterations the average error rate is: %f',\ %(numTests,errorRateSum/float(numTests)))
