Logistic回归

一、 预备知识

在学习之前，首先理解几个问题，有助于更好的完成接下来的阅读和python代码理解。

• 什么是回归？

  答：回归就是一种拟合过程。所谓拟合，就是指已知某函数的若干离散函数值$\{f_1,f_2,...f_n\}$，通过调整该函数中若干待定系数$f(\lambda_1,\lambda_2,...,\lambda_n)$，使得该函数与已知点集进行匹配，当达到最佳匹配时完成回归。

• 什么时候才能达到最佳匹配？

  答：该函数与已知点集的差别(最小二乘意义)最小时，达到最佳匹配。

• 有没有什么例外情况？

  答：这个可能就会涉及到过拟合、欠拟合以及剔除异常的情况，下面依次说明。

  • 过拟合：过拟合用字面意思就可以解释，即过度拟合。也就是说在拟合过程中为了完美匹配到所有的点而引入过于复杂的系数或者参数，和成语中的“矫枉过正”一个道理。

  • 欠拟合：欠拟合就是指没有达到应该达到的拟合程度，可能是由于模型建立的过于简单，或者参数引入过少导致。举个最简单的例子，散点具有二次函数的性质，但是做拟合时却只采用一次函数的模型进行拟合，这就是欠拟合。

  • 剔除异常：剔除异常，即在原始数据集中存在一些极度偏离总体数据的点，这些点称为异常点，而我们在拟合之前，进行数据预处理时，需要将这些点剔除。

• 有什么实用价值？

  答：拟合的最大价值就体现在分类上。

  闲扯：你进入一家专卖店，这时候营业员根据你身上的几个特征：上衣、裤子、鞋、颜色搭配、使用物品、说话方式等，就可以对你进行人为的分类。例如将你划分为学生、白领、程序员，诸如此类。当然，由于当今社会大数据的普及，已经不再需要营业员去人为的为你划分类别，而是根据摄像头和计算机，就能在一秒钟之内准确划分，并向你有针对性的推荐商品。

二、Logistic回归和Sigmoid函数的分类

1. 二值型输入分类器的数学原理

对于二值型输入分类器而言，最理想的情况是我们的输入数据经过函数处理后，能够划分为0-1两类，即处理函数的输入为0或1。

那么，有哪类函数可以完成这个目的呢？我们首先想到的一定是单位阶跃函数，在信号与系统中，记为$\epsilon$，具体如下图所示：

该函数具有十分优雅的性质，即将整个实数域的自变量映射到$[-1,1]$区间上。但是，我们会注意到，单位阶跃函数并不是连续函数，因此我们无法使用统一而简洁的函数表达式对其进行表示，这会让我们很不舒服。

闲扯：如果一个理论具有最优美的形式，那么它一定是具有最少附加条件的。这个假设和奥卡姆剃刀原理类似。另外，这个理论不仅仅适用于数学上，更在物理学上得到了最有力的验证。最有名的例子就是爱因斯坦的光电方程和麦克斯韦的电磁场理论，形式优美而简洁。爱因斯坦和薛定谔晚年试图统一场论，将全部理论用一个统一的方程进行阐述，也是对这种优雅理论的不懈追求，虽然最后未能成功。

这时候，就需要Sigmoid函数登场了，其具体计算公式如下：

$$\sigma(z)=\frac{1}{1+e^{-z}}$$

从图中可以看出，Sigmoid函数的映射就是将较大的输入值变得更大，较小的输入值变得更小，让两者区分得更加明显（注：此处大小不是指数值意义上的大小，而是指分类的明显程度），从而起到一个放缩的作用。因此，为了对数据进行分类，我们在使用Sigmoid函数时应该尽量避免出现输出值在0.5的较小邻域内出现，这样会影响我们最终的分类效果，好的输入应该尽量分布在靠近0或者靠近1的范围。

因此，为了实现Logistic回归分类器，我们可以在每个特征熵都乘以一个回归系数，然后把所有结果值相加，将这个综合代入Sigmoid函数中，进而得到一个0~1范围之间的数值。任何大于0.5的输入被分入1类，小于0.5的数据被分入0类。

所以，现在的问题可以转化为，最佳回归系数应该怎么求？

2. 如何求解最佳回归系数？

我们将Sigmoid函数的输入记为 $Z$，则输入和输出的关系用如下线性方程表示：（注：在这里我们不考虑输入和输出呈非线性的关系）

$$Z=w_0x_0+w_1x_1+w_2x_2+...+w_nx_n$$
采用向量形式，上述公式可以表示为：

$$Z=W^TX$$
其中， $Z,W,X$都是长度为 $n$的向量。

因此，我们所求解的最佳回归系数就是向量$W$的值，为了求解该向量，需要用到 梯度上升法。

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

那么，要沿着梯度方向进行寻找，我们需要去设定我们每次前进的距离，这就是步长$\alpha$。所以，用向量形式表示梯度上升算法中回归系数的更新为：

$$W=W+\alpha\nabla_wf(W)$$
该公式将一直迭代执行，知道达到某个停止条件为止。

3. 实现代码

说明：代码参考自《机器学习实战》，并针对python3进行了相关调整。

1. 回归参数求解

from numpy import *
import numpy as np

# ===============参数训练================
def loadDataSet():
    import re
    # 定义数据集和标签集
    dataMat = []; labelMat = []
    fr = open('testSet.txt')
    for line in fr.readlines():
        # 去掉空格，并将每一行的数据提取出来
        # 写成了fr.readline()，只能读取一个字符，为此debug了很久。被自己蠢哭。
        lineArr = line.strip().split()
        # print(lineArr)
        # 下行代码需要详细解释，见备注1-A
        dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])
        # 使用正则表达式匹配字符串，此处代码有误
        '''tempList1 = []; tempList2 = []
        tempList1 = lineArr[0]
        tempList2 = lineArr[1]
        print(tempList1)
        print(tempList2)
        if(tempList1[0] == '-'):
            print(re.match(r'\d*',tempList1).group(0))
            # num1 = float(0.0 - float(re.match(r'\d*',lineArr[0]).group(0)))
        else:
            num1 = float(lineArr[0])
        if(tempList2[0] == '-'):
            print(re.match(r'\d*', tempList2).group(0))
            # num2 = float(0.0 - float(re.match(r'\d*',lineArr[1]).group(0)))
        else:
            num2 = float(lineArr[1])
        dataMat.append([1.0, num1, num2])'''
        # 提取标签信息
        labelMat.append(int(lineArr[2]))
    return dataMat, labelMat

def sigmoid(inX):
    # 将结果值代入，计算sigmoid函数的值
    return 1.0/(1+exp(-inX))

# 梯度上升算法函数（重要指数：*****）
def gradAscent(dataMatIn, classLabels):
    dataMatrix = mat(dataMatIn)
    # 将标签矩阵转置(原：1xN的行向量；现：Nx1的列向量)
    labelMat = mat(classLabels).transpose()
    # 计算矩阵的行(m)和列(n)
    m, n = shape(dataMatrix)
    # 设定步长，alpha
    alpha = 0.001
    # 最大迭代次数为500
    maxCycles = 500
    # 初始化系数向量为元素为1的列向量
    weights = ones((n,1))
    for k in range(maxCycles):
        # h为nx1的列向量，用来存放分类得到的标签值
        h = sigmoid(dataMatrix*weights)
        # 正确标签和计算得到的标签之间的误差
        error = (labelMat - h)
        # 梯度上升，注意此处需要矩阵转置，是因为矩阵乘法需要行列对应。
        weights = weights + alpha * dataMatrix.transpose() * error
        # 以上三行公式的证明，见备注1-B

    return weights

2. 原始梯度上升算法的图形绘制

# =============图像绘制============
def plotBestFit(weights):
    # 使用python3画图时需要添加以下两行语句，否则报错。
    import matplotlib
    matplotlib.use('Qt5Agg')

    import matplotlib.pyplot as plt

    dataMat, labelMat = loadDataSet()
    # 将dataMat转化为数组
    dataArr = array(dataMat)
    # 取数组的行数作为数据个数
    n = shape(dataArr)[0]
    # 标签为1的数据的x、y坐标和标签为0的数据的x、y坐标
    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()
    # 111表示画一个大图
    ax = fig.add_subplot(111)
    # 将标签为1的点绘制为数量为30(s=30)，颜色为红色(c=red)，形状为正方形(marker='s')
    # 更多scatter用法：http://blog.youkuaiyun.com/u013634684/article/details/49646311
    ax.scatter(xcord1, ycord1, s = 30, c = 'red', marker = 's')
    ax.scatter(xcord2, ycord2, s = 30, c = 'green', marker = 'o')
    x = arange(-3.0,3.0,0.1)
    # 纵坐标的取值，具体解释见备注2-A
    y = (-weights[0]-weights[1]*x)/weights[2]
    ax.plot(x,y)
    plt.xlabel('X1'); plt.ylabel('X2')
    plt.show()
# =========测试程序========
dataArr, labelMat = loadDataSet()
weights = gradAscent(dataArr, labelMat)
# 注：此处使用getA()是为了返回一个n维数组，因为在计算y时用到了weights[1]*x，如果weights[1]不是一个n维数组，而是一个数时，维度不匹配。
plotBestFit(weights.getA())

3. 第一次改进的梯度上升算法

# =========算法升级=========
# 原算法在每次更新回归系数时需要遍历整个数据集，当样本量很大或者特征量很大时，计算复杂度过高！
# 改进方法是一次仅用一个样本点来更新回归系数，该方法称为随机梯度上升算法。
# 采用随机梯度上升算法对分类器进行增量式更新。
# 该方法是一个在线学习算法，与之对应的原梯度上升算法一次性处理全部数据，为批处理方式

def stocGradAscent0(dataMatrix, classLabels):
    m, n = shape(dataMatrix)
    alpha = 0.01
    weights = ones(n)
    # 循环部分和原梯度上升算法不同：
    # 1. 循环次数不再是设定的迭代次数-maxCycles，而是dataMatrix矩阵的行数，即样本个数
    for i in range(m):
        # 2. 利用回归方程计算结果值时，所求得结果为【数值】而非【向量】！
        # 注意：假定dataMatrix[i]包含两个特征x,y（即n=2），则梯度上升算法计算得到
        #     【（x*w1，y*w2）】,而随机梯度上升算法计算得到【x*w1+y*w2】
        h = sigmoid(sum(dataMatrix[i]*weights))
        # 3. error计算也是数值，而非向量。理由同上。
        error = classLabels[i] - h
        weights = weights + alpha * error * dataMatrix[i]
    return weights

# =========测试算法==========
# 测试算法得到的图可以说明，该随机梯度上升算法分类准确性下降明显，分类器错了接近1/3的样本。
dataArr, labelMat = loadDataSet()
weights = stocGradAscent0(array(dataArr), labelMat)
plotBestFit(weights)

4. 第二次改进的梯度上升算法

# =========算法再次升级=========
# 1. 改进后算法考虑了迭代次数的影响
# 2. 改进后的算法考虑了alpha值的调整，随迭代次数减小，降低了高次迭代时数据的波动
# 3. 改进后的算法在降低alpha的函数中，alpha每次减少1/(i+j)，其中j是迭代次数，i是样本点下标。
#    当j<<max(i)时，alpha就不是严格下降的。
# 其中第3条的优点需要补充：
def stocGradAscent1(dataMatrix, classLabels, numIter = 150):
    m, n = shape(dataMatrix)
    weights = ones(n)
    for j in range(numIter):
        dataIndex = list(range(m))
        for i in range(m):
            alpha = 4/(1.0+j+i) + 0.01
            # 随机选取样本点
            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

# ============测试算法============
weights = stocGradAscent1(array(dataArr), labelMat, 200)
plotBestFit(weights)

5. 相关备注

备注1-A：代码dataMat.append([1.0, float(lineArr[0]), float(lineArr[1])])中，append()包含了三项，其中后两项是数据的x坐标和y坐标，而第一项1.0用来构成常数项，即最终表达式为：

$$Z=w_0+w_1x_1+w_2x_2$$

备注1-B：

代码

 h = sigmoid(dataMatrix*weights)
 # 正确标签和计算得到的标签之间的误差
 error = (labelMat - h)
 # 梯度上升，注意此处需要矩阵转置，是因为矩阵乘法需要行列对应。
 weights = weights + alpha * dataMatrix.transpose() * error

为什么dataMatrix.transpose() * error可以表示梯度值？

解释如下：因为$Z=W^TX$，为了评估回归参数$W$的估计性能，我们定义Loss Function:

$$J(w)=\frac{1}{2}\Sigma(h_w(x^{(i)})-y^{(i)})^2$$
对上式求导（取梯度）得到：

$$\frac{dJ}{dw}=h_w(x)-y$$
其中， $h_w(x)$对应于代码中的 h = sigmoid(dataMatrix*weights)； $y$对应于代码中的labelMat。

备注2-A：

y = (-weights[0]-weights[1]*x)/weights[2]

因为$0=w_0+w_1x+w_2y$，调整 $y$的位置，有$y=(-w_0+w_1x)/w_2$。其中， $x,y$均为向量。