机器学习算法四——Logistic回归（1）

本章内容
⚫Sigmoid函数和Logistic回归分类器
⚫最优化理论初步
⚫梯度下降最优化算法
⚫数据中的缺失项处理

假设现在有一些数据点，我们用一条直线对这些点进行拟合（该线称为最佳拟合直线），这个拟合过程就称为回归。

Logistic回归进行分类的主要思想：根据现有数据对分类边界线建立回归公式，以此进行分类。

训练分类器的做法就是寻找最佳拟合参数，使用的是最优化算法。

一、基于Logistic回归和Sigmoid函数的分类
Logistic回归：计算代价不高，易于理解和实现；容易欠拟合，分类精度可能不高。
我们想要的函数是：接受所有的输入，然后预测出类别。例，在两个类时，函数输出0或1。——单位阶跃函数
单位阶跃函数的问题：该函数在跳跃点上从0瞬间跳跃到1，这个瞬间跳跃过程有时很难处理。
幸好，另一个函数也有类似的性质，就是Sigmoid函数：
σ(z) =1/(1+e^–z) .
下图给出了Sigmoid函数在不同坐标尺度下的两条曲线图：

如果横坐标刻度足够大，Sigmoid函数看起来很像一个阶跃函数。
为了实现Logistic回归分类器，我们可以在每个特征上都乘以一个回归系数，然后把所有的结果值相加，将这个总和代入Sigmoid函数中，进而得到一个范围在0~1之间的数值。所以，Logistic回归也可以被看成是一种概率估计。
所以，最佳回归系数是多少？如何确定它们的大小？

二、基于最优化方法的最佳回归系数确定
Sigmoid函数的输入记为z：
z = w₀x₀ + w₁x₁ + w₂x₂ +……+ w_nx_n
向量写法： z=w^Tx
向量x是分类器的输入数据，向量w是我们要找到的最佳参数。为寻找最佳参数，下面介绍几种方法。
1、梯度上升法
思想：要找到某函数的最大值，最好的方法是沿着该函数的梯度方向探寻。如果梯度记为▽，则函数f(x,y)的梯度由下式表示：

这个梯度意味着要沿x的方向移动，沿y的方向移动。其中，函数f(x,y)必须要在待计算的点上有定义并且可微。

除了移动方向，还要考虑移动量的大小，该量值称为步长，记作α。梯度算法的向量迭代公式：
该公式将一直被迭代执行，直至达到某个停止条件为止，比如迭代次数达到某个指定值或算法达到某个可以允许的误差范围。
梯度上升算法用来求函数的最大值，而梯度下降算法用来求函数的最小值：

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

上图中有100个样本点，每个点包含两个数值型特征：X1和X2。在此数据集上，我们将通过使用梯度上升法找到最佳回归系数，也就是拟合出Logistic回归模型的最佳参数。
下面是梯度上升算法的具体实现：

import numpy as np

#Logistic回归梯度上升优化算法
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])])#将每一行的两个特征x1，x2，加上x0=1,组成列表并添加到数据集列表中
        labelMat.append(int(lineArr[2]))#将当前行标签添加到标签列表
    return dataMat,labelMat #返回数据列表，标签列表

#定义sigmoid函数
def sigmoid(inx):
    return 1.0/(1+ np.exp(-inx))

#梯度上升法更新最优拟合参数
def gradAscent(dataMatIn,classLabels):#@dataMatIn：数据集  @classLabels：数据标签
    dataMatrix=np.mat(dataMatIn) #将数据集列表转为Numpy矩阵
    labelMat=np.mat(classLabels).transpose() #将标签列表转为Numpy矩阵，并转置
    m,n=np.shape(dataMatrix) #获取数据集矩阵的行数和列数
    alpha=0.001 #学习步长
    maxCycles=500 #最大迭代次数
    weights=np.ones((n,1)) #初始化权值参数向量每个维度均为1.0
    
    for k in range(maxCycles): #循环迭代次数     
        h=sigmoid(dataMatrix*weights) #求当前的sigmoid函数预测概率
        error=(labelMat-h) #此处计算真实类别和预测类别的差值                                                   
        weights=weights+alpha*dataMatrix.transpose()*error #按照差值的方向调整回归系数

    return weights

def main():
    dataArr, labelMat = loadDataSet()
    weights = gradAscent(dataArr, labelMat)
    print(weights)

顺便查了下import numpy 和 from numpy import *的区别：尽量使用前者，后者容易出现命名上的混淆，前者使用时加上前缀（比如np.ones((n,1))）。
返回的是回归系数：

>>>[[ 4.12414349]
 [ 0.48007329]
 [-0.6168482 ]]

3、分析数据：画出决策边界
解出的回归系数确定了不同类别数据之间的分隔线，下面我们画出分隔线，使得优化的过程便于理解。
代码如下：

#画出数据集和Logistic回归最佳拟合直线的函数
def plotBestFit(wei):  
   import matplotlib.pyplot as plt  
   weights = wei.getA()   #getA()是numpy的一个函数
   dataMat, labelMat = loadDataSet()  
   dataArr = np.array(dataMat)  
   n = np.shape(dataArr)[0]  #数据集的行数，即样本数
   xcord1 = []; ycord1 = []  #一类数据的x、y坐标（x、y坐标代表两个特征值）
   xcord2 = []; ycord2 = []  #另一类数据的x、y坐标
   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)  #-3到3，步长为0.1
   y = (-weights[0]- weights[1]*x)/weights[2]  #最佳拟合直线
   ax.plot(x, y)  
   plt.xlabel('X1');  #横坐标代表X1
   plt.ylabel('X2');  #纵坐标代表X2
   plt.show()  

def main():
   dataArr, labelMat = loadDataSet()
   weights = gradAscent(dataArr, labelMat)
   plotBestFit(weights)

说明一下这个最佳拟合直线
y = (-weights[0]- weights[1]*x)/weights[2]的由来：
由于0是两个分类（类别0和类别1）的分界处，因此，设定0= w₀x₀ + w₁x₁ + w₂x₂ ，解出x1和x2的关系式（即分隔线的方程）：
x₂ = （-w₀-w₁x₁）/ w₂

结果如下：

这个分类结果相当不错，从图上看只错分了两到四个点。但是，尽管例子简单且数据集很小，此方法却需要大量的计算（300次乘法）。
下面，对该算法稍作改进。

4、训练算法：随机梯度上升
梯度上升算法每次更新回归系数时，都需要整个数据集。该方法处理100个左右的数据集时尚可，但如果有数十亿样本和成千上万的特征，那计算的复杂度实在太高。
一种改进方法是一次仅用一个样本点来更新回归系数，称为随机梯度上升算法。
由于可以在新样本到来时对分类器进行增量式更新，因而随机梯度上升算法是一个在线学习算法 。与在线学习相对应，一次处理所有数据被称作是==“批处理“==。
下面是代码实现：

#随机梯度上升算法
def stocGradAscent0(dataMatrix, classLabels):
    m,n = np.shape(dataMatrix)
    alpha = 0.01
    weights = np.ones(n)
    for i in range(m):
        h = sigmoid(sum(dataMatrix[i]*weights))
        error = classLabels[i] - h
        dataMatrix[i]=np.mat(dataMatrix[i])
        weights = weights + alpha * error * dataMatrix[i]   
    return weights

def main():
    dataArr, labelMat = loadDataSet()
    weights = stocGradAscent0(np.array(dataArr), labelMat)
    plotBestFit(weights)

这里运行后会报错：

>>>weights = wei.getA()  #getA()是numpy的一个函数
AttributeError: 'numpy.ndarray' object has no attribute 'getA'

是因为在输入时已经转换为数组了，这里需要注释掉plotBestFit()的矩阵转为数组：

#weights = wei.getA()  #getA()是numpy的一个函数
 weights = wei

最终：

可以看到拟合出来的直线不像上一张图那么完美，错分了1/3的样本。
直接比较代码并不公平，后者的结果是在数据集上迭代了500次得到的。一个判断优化算法优劣的可靠方法是看它是否收敛，即参数是否达到了稳定值，是否还会不断地变化？
下面改进随机梯度上升算法，使其在整个数据集上运行200次。代码如下：

#改进的随机梯度上升算法
def stocGradAscent1(dataMatrix, classLabels, numIter=150):
    m,n = np.shape(dataMatrix)
    weights = np.ones(n)   #initialize to all ones
    for j in range(numIter):
        dataIndex = list(range(m))
        for i in range(m):
            alpha = 4/(1.0+j+i)+0.01   #alpha每次迭代时需要调整
            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

def main():
    dataArr, labelMat = loadDataSet()
    weights = stocGradAscent1(np.array(dataArr), labelMat, numIter=500)
    plotBestFit(weights)

该进一：alpha在每次迭代的时候都会调整，虽然alpha会随着迭代次数不断减小，但永远不会减小到0，因为其中还存在一个常数项。（保证在多次迭代之后新数据仍有一定的影响）
若处理的问题是动态变化的，可适当加大上述常数项，确保新的值获得更大的回归系数。alpha每次减少1/(j+1)，当j<<max(i)时，alpha就不是严格下降的。
改进二：通过随机选取样本来更新回归系数。这种方法将减少周期性的波动。

从结果可以看到，分隔线达到了与GradientAscent()差不多的效果，但是所使用的计算量更少。

下面我们需要完成具体的分类任务，使用随机梯度上升算法来解决病马的生死预测问题。