Logistic回归算法

前言

关于Logistic的重要性，不言而喻，网上也有许多关于Logistic的资料，这篇文章我就写写我自己在看了机器学习实战后，对于Logistic的理解吧。首先我们要知道Logistic能干什么？作为我们学习的第一个回归算法，它根据现有的数据对分类边界建立回归公式，以此进行分类。我们的目的就是找到那根最佳拟合直线。也就是我们能把数据分成两类，那么我们便可以把这应用在一些数值的分类问题上，比如怎么判断图片里内容是不是猫（如下图），怎么预测一匹马是否能被治愈等二值问题，以及后面我们深度学习神经网络里面就是以Logistic为神经基础单元，用来图像识别，语言识别等等。
可能你会发现与k-近邻算法有点像，但是k-近邻算法能有不止两类的分类，它是根据离分类数据点的距离多数表决来进行分类。

算法描述

这里我们要学习一个Sigmoid函数，计算公式如下：

σ(z)=11+e−z

我们这里可以根据matplotlib画出此函数的图。

import sysfrom pylab import *import matplotlib.pyplot as pltimport numpy as npt = np.arange(-60.0,60.3,0.1)s = 1/(1 + np.exp(-t))ax = plt.subplot(211)ax.plot(t,s)ax.axis([-5,5,0,1])plt.xlabel('x')plt.ylabel('Sigmoid(x)')ax = plt.subplot(212)ax.plot(t,s)ax.axis([-60,60,0,1])plt.xlabel('x')plt.ylabel('Sigmoid(x)')plt.show()

我们调整横坐标为不同的尺度可以看到此函数在[0，1]范围内，正好符合概率值范围，这个函数有利于我们算法计算（这只是算法的一种函数选择，我们也可以选择其它函数，但是效果不见得那么好），二值即为0,1分类，我们需要预测为0,1的概率。
对于任意x∈R，y∆ = p（y=1|x），y∆∈[0,1]，此式代表的是任意给一个x，y=1的概率。
由于是线性回归，我们联立三个公式：

⎧⎩⎨⎪⎪⎪⎪⎪⎪z=w0x0+w1x1+w2x2+...+wnxnσ(z)=11+e−zerror=y−y′

y代表数据实际分类值，y’代表预测分类值概率，也就是σ（z）。
上面三个式子我们如果用高中数学肯定解不开，未知量过多，这里我们可以采用优化的方法来使得error最小，从而求出最佳参数w拟合值。

• 基于梯度上升法的优化方法
  其实梯度上升法与我们经常讲的梯度下降法差不多，只不过前者求的是函数的最大值，后者求的是函数最小值，公式中符号不同而已。公式为：
  
  w:=w+α∆f(w)
  
  怎么理解这个公式呢？画个图就明白了。
  
  当w>0时，可以看到导数小于0，公式中∆f(w)<0，α控制梯度的步长>0,：=是相对于程序中的赋值语句，所以w会减小，值会向函数最大值方向移动。
  当w<0时，可以看到导数大于0，公式中∆f(w)>0，α控制梯度的步长>0,：=是相对于程序中的赋值语句，所以w会增大，值会向函数最大值方向移动。
  最后函数会无限接近最大值，达到某个可允许误差范围，直至收敛。
  如果是二维的，我们一样按照偏导的方向去接近最大值。
  

训练算法：使用梯度上升找到最佳参数

• 这里我们用个例子来训练logistic算法，并用梯度上升来找到最佳拟合参数。
  数据集

import numpy as npdef loadDataSet():#导入数据集    dataMat = []; labelMat = []    fr = open('C:\\Users\\user\\Desktop\\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,labelMatdef sigmoid(inX):#sigmoid函数实现    return 1.0/(1+np.exp(-inX))def gradAscent(dataMatIn, classLabels):#梯度上升法选取最优参数    dataMatrix = np.mat(dataMatIn)                 labelMat = np.mat(classLabels).transpose()     m,n = np.shape(dataMatrix)    alpha = 0.001    maxCycles = 500    weights = np.ones((n,1))    for k in range(maxCycles):                      h = sigmoid(dataMatrix*weights)             error = (labelMat - h)                      weights = weights + alpha * dataMatrix.transpose()* error #matrix mult    return weights

数据集是100行3列，最后一列为数据的类别，这里我们先处理数据集，将数据集第一二列导入列表，我们还加入一列全为1.0的列，作为x0的输入。
第三个函数，我们通过误差的正负来修正weights。

这里输出迭代500次的weights的值。

• 分析数据：画出决策边界

def plotBestFit(weights):    import matplotlib.pyplot as plt    dataMat,labelMat=loadDataSet()    dataArr = np.array(dataMat)    n = np.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 = np.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()

上面就是把数据坐标放入列表里面，然后通过scatter（）函数画散点图。

 plotBestFit(weights.getA())

这里画的时候，我们要把矩阵转化为数组，因为plot（）函数只接受数组。

这里我们看到，对于分类结果还是比较满意的（除了一个异常点）。

• 训练算法：随机梯度上升
  随机梯度上升就是解决当有数十亿样本和成千上万的特征，那么该方法的计算复杂度就有点高了。随机梯度上升一次仅用一个样本点来更新回归系数，这有什么好处呢？这可以方便我们实时更新回归系数，算是一个“在线学习”算法。
  这里我们修改stocGradAscent0函数

def stocGradAscent1(dataMatrix,classLabels,numIter=150):    m,n = np.shape(dataMatrix)    weights = np.ones(n)    for j in range(numIter):        dataIndex = list(range(m))        for i in range(m):            alpha = 4/(1.0+i+j)+0.01            randIndex = int(np.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步长也随着迭代次数在减小，方便最后收敛的值更精确，这会缓解最后的数据波动。
随机选取样本可以缓解周期性波动。
我们先来看看改正后的函数的收敛图：

可以看到收敛速度特别快，大约在5次遍历所有数据集时，就已经收敛了。
我们再看看回归拟合的怎么样：

weights = stocGradAscent1(array(dataArr),labelMat,500)plotBestFit(weights)

可以看到并不比批处理差，何况速度还更快，何乐而不为呢。

示例：从疝（shan）气病症预测病马的死亡率

这个数据集中有30%的值是缺失的，一般只有两种方法可选：

• 删除含有缺失数据的实例
• 用合理的替代值替换缺失值（均值、-1、0，用相似样本的均值、预测缺失值）

这里我们就用0来代替所有缺失值，这样在更新w时也不会产生影响。

1. 测试算法

def classifyVector(inX,weights):    prob = sigmoid(sum(inX*weights))    if prob > 0.5:return 1.0    else:return 0.0def colicTest():    frTrain = open('C:\\Users\\user\\Desktop\\horseColicTraining.txt')    frTest = open('C:\\Users\\user\\Desktop\\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(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(float(currLine[i]))        if int(classifyVector(np.array(lineArr),trainWeights))!=int(currLine[21]):            errorCount += 1    errorRate = (float(errorCount)/numTestVec)    print("错误率为：%f" % errorRate)    return errorRatedef multiTest():    numTests = 10;errorSum = 0.0    for k in range(numTests):        errorSum += colicTest()    print("10次平均后的错误率为:%f" %(errorSum/float(numTests)))multiTest()

我们首先创建了一个分类函数，如果概率大于0.5的，我们分类为1，小于则为0。
然后采用随机梯度上升，迭代500次，最后得出错误率。我们取十次误差的平均值得出结果：

这是我的结果，你的结果稍有不同，不过我认为错误率35%有点高了，与实际应用还是有差距的。