import matplotlib.pyplot as plt from numpy import * def loadData(): dataMat=[] labelMat=[] fr = open('sample1.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])) # print(dataMat) # print(labelMat) return dataMat, labelMat def plotDot(): dataMat,labelMat=loadData() dataArr = array(dataMat) n = shape(dataArr)[0] x0=[] x1=[] y0=[] y1=[] for i in range(n): if labelMat[i]==1: x1.append(dataMat[i][1]) y1.append(dataMat[i][2]) else: x0.append(dataMat[i][1]) y0.append(dataMat[i][2]) print(x0) print(x1) fig = plt.figure() ax = fig.add_subplot(111)#分割为1行1列第1块 ax.scatter(x0, y0, s=30, c='red', marker='s') ax.scatter(x1, y1, s=30, c='green') plt.show() plotDot()
数据集中一共有100个点，每个点包含两个数值型特征：X1和X2。因此可以将数据在一个二维平面上展示出来。我们可以将第一列数据(X1)看作x轴上的值，第二列数据(X2)看作y轴上的值。而最后一列数据即为分类标签。根据标签的不同，对这些点进行分类。
梯度上升法的伪代码：
每个回归系数初始化为1 重复R次： 计算整个数据集的梯度 使用alpha * gradient更新回归系数的向量 返回回归系数
代码实现：
from numpy import * #预处理数据 def loadDataSet(): #创建两个列表 dataMat=[];labelMat=[] #打开文本数据集 fr=open('sample1.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(dataMatrix) #学习步长 alpha=0.001 #最大迭代次数 maxCycles=500 #初始化权值参数向量每个维度均为1.0 weights=ones((n,1)) print(weights) #循环迭代次数 for k in range(maxCycles): #求当前的sigmoid函数预测概率 h=sigmoid(dataMatrix*weights) #*********************************************** #此处计算真实类别和预测类别的差值 #对logistic回归函数的对数释然函数的参数项求偏导 error=(labelMat-h) #更新权值参数 weights=weights+alpha*dataMatrix.transpose()*error #*********************************************** return weights
dataMat,labelMat=loadDataSet() m=gradAscent(dataMat,labelMat) print(m)
我们知道对回归系数进行更新的公式为：
其中
是对参数w求偏导数。则我们可以通过求导验证logistic回归函数对参数w的梯度为(y_{i}-sigmoid
5.2、画出决策边界
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() plotBestFit(m)
在拟合分类线中，我们设定sigmoid函数为0。0是两个分类的分解处，因此我们设定
0=w_{0}x_{0}+w_{1}x_{1}+w_{2}x_{2},然后根据函数制定x坐标表示为x1，y坐标表示为x2，画出分类线。尽管分类结果不错，但是这个方法却需要大量计算。
5.3、 随机梯度上升
梯度上升算法在每次更新回归系数时需要遍历整个数据集，该方法在处理100个左右的数据集时尚可，但如果有数十亿样本和成千上万的特征，那么该方法就很费劲了。所有我们思考了优化，一次仅用一个样本点来更新回归系数，该方法称为随机梯度上升法。对应的更新公式是：
由于可以在新的样本到来时对分类器进行增量式更新，因而随机梯度上升算法是一个在线学习算法。与之对应的是批处理算法。
随机梯度上升伪代码：
所有回归系数初始化为1 对数据集中每个样本 计算该样本的梯度 使用alpha*gradient更新回归系数值 返回回归系数值
代码：
def stocGradAscent0(dataMatrix,classLabels): m,n=shape(dataMatrix) alpha=0.01 weights=ones(n) dataMatrix=array(dataMat) for i in range(m): h=sigmoid(sum([dataMatrix[i]*weights])) error = classLabels[i]-h weights=weights+alpha*error*dataMatrix[i] print(type(dataMatrix[i])) return weights
随机梯度上升法和批量梯度上升法是两个极端，一个采用所有数据来梯度上升，一个用一个样本来梯度上升。自然各自的优缺点都非常突出。对于训练速度来说，随机梯度上升法由于每次仅仅采用一个样本来迭代，训练速度很快，而批量梯度上升法在样本量很大的时候，训练速度不能让人满意。对于准确度来说，随机梯度上升法用于仅仅用一个样本决定梯度方向，导致解很有可能不是最优。对于收敛速度来说，由于随机梯度上升法一次迭代一个样本，导致迭代方向变化很大，不能很快的收敛到局部最优解。 但值得一提的是，随机梯度上升法在处理非凸函数优化的过程当中有非常好的表现，由于其上升方向具有一定随机性，因此能很好的绕开局部最优解，从而逼近全局最优解。
一个判断优化算法的可靠办法是看它是否收敛，也就是说参数是否达到了稳定值，是否还会变化。对此观察随机梯度上升算法的3个参数的变化：
def stocGradAscent0(dataMatrix,classLabels): m,n=shape(dataMatrix) print(n) x1 = arange(0, 142) y1=[] y2=[] y3=[] alpha=0.01 weights=ones(n) dataMatrix=array(dataMat) for i in range(142): h=sigmoid(sum([dataMatrix[i]*weights])) error = classLabels[i]-h weights=weights+alpha*error*dataMatrix[i] y1.append(weights[0]) y2.append(weights[1]) y3.append(weights[2]) fig=plt.figure(figsize=(14,14),dpi=100) plt.subplot(3, 1, 1) plt.plot(x1,y1) plt.subplot(3, 1, 2) plt.plot(x1, y2) plt.subplot(3, 1, 3) plt.plot(x1, y3) # ax = fig.add_subplot(111) # bx = fig.add_subplot(211) # cx = fig.add_subplot(311) # ax.plot(x1, y1) # bx.plot(x1,y2) # cx.plot(x1,y3) plt.show() return weights
我们发现第三个参数仅仅在50次迭代就达到了稳定值，而其他两个则需要更多次迭代。产生这种现象的原因是存在一些不能正确分类的样本点（数据集并非线性可分），在每次迭代时会引发系数的剧烈变化。
5.4、优化随机梯度上升
我们对随机梯度上升算法做一些调整：
def stocGradAscent0(dataMatrix,classLabels,numIter=150): m,n=shape(dataMatrix) alpha=0.01 weights=ones(n) dataMatrix=array(dataMat) for j in range(150): dataIndex = range(m) dataIndex=list(dataIndex) for i in range (m): #alpha每次迭代时调整 alpha=4/(1.0+j+i)+0.1 #随机选取更新 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
该分割线达到了与gradAscent（）差不多的效果，而且计算量明显少于批量梯度上升算法。
六、实战：从疝气病症预测病马的死亡率
我们将使用Logistic回归来预测患有疝气病的马存活问题。此数据集包含368个样本和28个特征。另外需要说明的是这些数据集是有部分缺失，我们需要先进行数据预处理再利用回归预测。
对于数据集中缺失数据我们可以采取以下办法进行处理：
- 使用可用特征的均值来填补缺失值；
- 使用特殊值来填补缺失值，比如-1；
- 忽略有缺失值的样本；
- 使用相似样本的均值添补缺失值；
- ‘使用另外的机器学习算法预测缺失值；（如KNN）
而对于该数据集我们选择用0代替缺失值，因为根据回归系数的更新公式：
当一数据集特征为0时，
将不作更新。
另外由于sigmoid（0）=0.5，它对结果的预测不具有任何倾向性，因此上述做法也不会对误差项造成任何影响。
现在我们根据上述所做的工作，现在只需要将特征向量乘以最优化系数得到的值全部相加，作为输入给sigmoid函数中。如果sigmoid值大于0.5则给定标签为1，否则为0。
def colicTest(): frTrain = open('horseColicTraining.txt', encoding='utf8'); frTest = open('horseColicTest.txt', encoding='utf8') trainingSet = [] trainingLabels = [] print(frTest) print(frTrain) 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 = stocGradAscent0(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("该模型错误率为: %f" % errorRate) return errorRate
由于每次是随机取数，所以我们可将colicTest()迭代次数更多一些：
def multiTest(): numTests = 10 errorSum = 0.0 for k in range(numTests): errorSum += colicTest() print("该模型迭代%d次平均错误率为: %f",numTests,(numTests, errorSum / float(numTests)))
运行输出：
可见错误率还是很低，因为数据存在30%的丢失。
总结
Logistic回归利用Sigmoid函数作为最佳拟合函数，利用最优化算法求出参数，再把数据集特征向量矩阵与对应训练好的参数相乘再相加，作为Sigmoid的输入从而得到预测。
