机器学习第五章 logistic回归

引言

利用logistic回归进行分类的主要思想是：根据现有数据对分类边界线建立回归公式，以此进行分类。
Logistic 回归优缺点
优点：计算代价不高，易于理解和实现。
缺点：容易欠拟合，分类精度可能不高。
适用数据类型：数值型和标称型数据。
Logistic 回归的一般过程
1 收集数据：采用任意方法收集数据。
2 准备数据：由于需要进行距离计算，因此要求数据类型为数值型。另外，结构化数据格式则最佳。
3 分析数据：采用任意方法对数据进行分析。
4 训练算法：大部分时间将用于训练，训练的目的是为了找到最佳的分类回归系数。
5 测试算法：一旦训练步驟完成，分类将会很快。
6 使用算法：首先，我们需要输入一些数据，并将其转换成对应的结构化数值；接着，基于训练好的回归系数就可以对这些数值进行简单的回归计算，判定它们属于哪个类别，在这之后，我们就可以在输出的类别上做一些其他分析工作。

5.1基于logistic回归和sigmoid函数的分类

我们想要的函数是，能接受所有的输入然后预测出类别。该函数称为海维赛德阶跃函数，或者直接称为单位阶跃函数。然而，该函数存在一些问题：在跳跃点上从0瞬间跳跃到1，这个瞬间跳跃过程有时很难处理。而另一个函数也有类似的性质，且在数学上更易处理，它就是Sigmoid函数。公式如下：

Sigmoid函数图如图，当横坐标刻度足够大时，Sigmoid函数看起来就像阶跃函数一样。

因此，为了实现Logistic回归分类器，我们可以在每个特征上都乘以一个回归系数，然后把所有的结果值相加，将总和带入Sigmoid函数中，得到一个范围在0～1之间的数值。大于0.5的数据被分入1类，小于0.5被归入0类。所以，Logistic回归也可以被看成是一种概率估计。

5.2基于最优化的最佳回归系数确定

Sigmoid函数的输入记为z ，由下面公式得出：
$$z=w0x0+w1x1+w2x2+...+wnxn$$
如果采用向量的写法，上述公式可以写成z=wT其中的向量x是分类器的输入数据，向量w也就是我们要找到的最佳参数（系数），从而使得分类器尽可能的精确。

5.2.1梯度上升法

梯度上升法基于的思想是：要找到某函数的最大值，最好的方法是沿着该函数的梯度方向探寻。

实际就是函数f ff的偏导数。
函数f(x,y)必须要在待计算的点上有定义并且可微。

​梯度上升算法到达每个点后都会重新估计移动的方向。从P0开始，计算完该点的梯度，函数就根据梯度移动到下一点P1。在P1点，梯度再次被重新计算，并沿新的梯度方向移动到P2。如此循环迭代，直到满足停止条件。迭代的过程中，梯度算子总是保证我们能选取到最佳的移动方向

梯度算法的迭代公式如下：
$$w=w+\alpha \nabla wf((w)$$
α称为步长,代表移动量的大小。该公式将一直被迭代执行，直至达到某个停止条件为止，比如迭代次数达到某个指定值或算法达到某个可以允许的误差范围。

5.2.2训练算法

梯度上升法的伪代码如下：
    每个回归系数初始化为1
    重复R次：
        计算整个数据集的梯度
        使用alpha*gradient更新回归系数的向量
        返回回归系数

​新建 logRegres.py 文件，输入代码：

from numpy import  *
 
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)             #将矩阵转换成NumPy矩阵
    labelMat = mat(classLabels).transpose() #转置矩阵
    m,n = shape(dataMatrix)
    alpha = 0.001#步长
    maxCycles = 500#迭代次数
    weights = ones((n,1))
    for k in range(maxCycles):              #矩阵运算
        h = sigmoid(dataMatrix*weights)     #多次矩阵运算
        error = (labelMat - h)              #vector subtraction
        weights = weights + alpha * dataMatrix.transpose()* error #matrix mult
    return weights

代码在开头提供了一个便利函数 loadDataSet()，它的主要功能是打开文本文件 testSet.txt 并逐行读取。每行前两个值分别是X1和X2,第三个值是数据对应的类别标签。此外，为了方便计算,该函数还将X0的值设为1.0。接下来的函数是 sigmoid()。
梯度上升算法的实际工作是在函数 gradAscent () 里完成的，该函数有两个参数。第一个参数是dataMathIn，它是一个2维NumPy数组，每列分别代表每个不同的特征，每行则代表每个训练样本。我们现在采用的是100个样本的简单数据集，它包含了两个特征X1和X2，再加上第0维特征X0，所以dataMathln里存放的将是100×3的矩阵。在处，我们获得输人数据并将它们转换成NumPy矩阵。第二个参数是类别标签，它是一个1×100的行向量。为了便于矩阵运算，需要将该行向量转换为列向量，做法是将原向量转置，再将它赋值给labelMat。接下来的代码是得到矩阵大小，再设置一些梯度上升算法所需的参数。
变量alpha是向目标移动的步长，maxCycles是迭代次数。在for循环迭代完成后，将返回训练好的回归系数。需要强调的是，所做的运算是矩阵运算。变量h不是一个数而是一个列向量，列向量的元素个数等于样本个数，这里是100。对应地，运算dataMatrix * weights代表的不止一次乘积计算，事实上该运算包含了300次的乘积。

测试代码如下所示

import logRegres
dataArr,labelMat = logRegres.loadDataSet()
print(logRegres.gradAscent(dataArr,labelMat))

输出结果

5.2.3分析数据：画出边界线

为logRegres.py添加代码

def plotBestFit(weights):
    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()

代码中设置了sigmoid函数为0。0是两个分类（类别1和类别0）的分界处。因此，设定 0=w0x0 + w1x1 + w2x2，然后解出X1和X2的关系式（即分割线的方程，注意X0=1）。

测试代码如下所示

import logRegres
dataArr,labelMat = logRegres.loadDataSet()
weights = logRegres.gradAscent(dataArr,labelMat)
logRegres.plotBestFit(weights.getA())

输出结果

5.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数组。
测试代码如下所示

import logRegres
from numpy import *
dataArr,labelMat = logRegres.loadDataSet()
weights = logRegres.stocGradAscent0(array(dataArr),labelMat)
logRegres.plotBestFit(weights)

输出结果如下所示

拟合直线效果尚可，但分类效果并不好。两者算法效果虽然有差别，但迭代次数也不相同，前者迭代次数多于后者。一个判断优化算法优劣的可靠方法是看它是否收敛，也就是说参数是否达到稳定值，是否还会不断变化。

对此 ，我们在程序随机梯度上升算法上做了些修改，使其在整个数据集上运行200次。
下图图为梯度上升算法和改进的梯度上升算法中回归系数与迭代次数的关系。可以看出，改进后的随机梯度上升算法随着迭代次数的增加回归系数与迭代次数关系更加接近阶跃函数曲线。

5.3示例：从气病症预测病马的死亡率

5.3.1准备数据

当数据缺失时可选的做法：

1 使用可用特征的均值来填补缺失值；
2 使用特殊值来填补缺失值；
3 忽略有缺失值的样本；
4 使用相似样本的均值填补缺失值；
5 使用另外的机器学习算法预测缺失值

我们对下一节要用的数据集进行预处理，使其可以顺利地使用分类算法。在预处理阶段 需要做两件事：
第一，所有的缺失值必须用一个实数值来替换，因为我 们使用的numpy数据类型不允许包含缺失值。这里选择实数0来替换所有缺失值，恰好能适用回归。

第二，如果在测试数据集中发现了一条数据的类别标签已经缺失，那么我们的简单做法是将该条数据丢弃。这是因为类别标签与特征不同，很难确定采用某个合适的值来替换。

5.3.2 测试算法：用Logistic回归进行分类

使用Logistic回归方法进行分类并不需要做很多工作，所需做的只是把测试集上每个特征向量乘以最优化方法得来的回归系数，再将成绩结果求和，最后输入到Sigmoid函数中即可。如果对应的Sigmoid值大于0.5就预测类别标签为1，否则为0。

将下面的代码添加到logRegres.py中

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) #n行22列
        trainingLabels.append(float(currLine[21]))
    trainWeights = stocGraAscent1(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'the error rate is {errorRate}')
    return errorRate

#多次测试，求错误率
def multiTest():
    numTest = 10
    errorSum= 0.0
    for k in range(numTest):
        errorSum += colicTest()
    print(f'after {numTest} iterations the average error rate is {errorSum/float(numTest)}')

5.4小结

本章引入了logistic回归，目的就是为了找到一个非线性函数SIgmoid的最佳拟合参数，求解过程可以由最优化算法来完成。在最优化算法中，最常用的就是梯度上升算法，而梯度上升算法又可以简化为随机梯度上升算法。此外引入了一些处理缺失数据的方法。