机器学习之Logistic回归

冒着被老板骂的风险，没有继续调试论文上的代码，学习了一个新的算法Logistic回归。之前已经在李航老师的《统计学习》上过了一遍原理，这次是具体的代码实现。

一、Logistic回归原理

机器学习算法分面临两大任务：分类与回归。分类就是将数据集划分成合适的类中，针对的是标称型数据，比如常见算法的有KNN算法、支持向量机（SVM）算法、决策分类树、朴素贝叶斯算法等。回归则是用于预测数值型数据，比如决策回归树。Logistic回归就是这么一种特殊的存在，看起来像是回归算法，其实是一种分类算法，只是用回归的思想解决了分类问题。

1.Sigmoid函数

之前说过，Logistic函数是用回归的思想来解决分类问题。既然有回归，就要有回归函数，理想的回归函数就是，我们给定输入，就会有相应的类别输出。在二项逻辑斯谛回归模型中，只有两类，分别为0,1。在这里，单位阶跃函数能解决这个问题，但直接这样设计阶跃函数不方便后面的优化计算，因为函数值不连续，无法进行一些相关求导。 所以，逻辑回归中，大家选了一个统一的函数，即Sigmoid函数，$$\sigma (z)=\frac{1}{1+e^{-z}}$$
主体思路是这样的：为了实现Logistic回归分类器，我们在每个特征上乘以一个回归系数，然后将所有的和相加，$$z=w_0{x}_0+w_1{x}_1+\cdots +w_n{x}_n$$
再将这个和代入到Sigmoid函数中，得到一个范围在0-1之间的数值。观察数值，将大于0.5的数据归为1类，小于0.5的归为0类。

则逻辑回归输出的预测函数数学表达式为 $$h_{\theta }(x)=g({\theta }^{\tau }x)=\frac{1}{1+e^{-{\theta }^{\tau }x}}$$

2.逻辑斯谛回归模型

这里，我们讨论的是二项逻辑斯谛回归模型。
定义：
$$P(Y=1|x)=\frac{exp(w\bullet x)}{1+exp(w\bullet x)}$$
$$P(Y=0|x)=\frac{1}{1+exp(w\bullet x)}$$
其中，$w$为权值向量，$x$为输入。
逻辑斯谛回归就是比较这两个条件概率的大小，将实例x分到概率值较大的一类。

3.逻辑回归和线性回归

我们知道，回归可以分为线性回归和逻辑回归，那么这两种方法有什么异同呢：
相同：

1. 都使用了极大似然估计来对模型参数进行估计。线性回归使用最小二乘法，实际上是在自变量$x$和参数$\theta$确定，因变量$y$服从正态分布的假设下，使用极大似然估计的一个化简。逻辑回归通过对极大似然估计的学习，学到最佳参数$\theta$
2. 在求解参数的过程中，都使用梯度下降或梯度上升算法。

不同：

1. 逻辑回归处理的是分类问题，线性回归处理的是回归问题。
2. 二项逻辑回归因变量的取值为一个二元分布，模型学习得到的是$E[y|x;\theta ]$;而线性规划实际上求解的是$y={\theta }^{\tau }x$,，是对假设的真实关系$y={\theta }^{\tau }x+ϵ$的一个近似。
3. 逻辑回归的因变量是离散的，而线性回归的因变量是连续的。

二、梯度上升算法和随机梯度上升算法

由Logistic原理可知，我们可以确定一个分类器的函数形式，那么问题来了：为了使分类器的效果达到最优，每一个特征上的回归系数应该为多少？ 使用梯度上升法

1.梯度上升法

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

从$x_0$开始，计算完该点的梯度，函数沿着得到的梯度移动到下一个点$x_1$，在$x_1$点梯度被重新计算，并沿着新的梯度方向移动到$x_2$，如此循环迭代，一直到满足迭代条件为止。
如果数据集比较小，完全可以采用全数据集（Full Batch Learning）的形式，采用全数据集有两个好处：

1. 由全数据集确定的方向能够更好的代表样本总体，从而更准确的朝着极值所在的方向移动；
2. 由于不同权重的梯度值差别巨大，因此选取一个全局的学习率很困难。Full Batch Learning 可以使用Rprop只基于梯度符号并且针对性单独更新各权值。

类比梯度下降法，梯度上升法存在一些缺点：

1. 梯度下降法并不能保证被优化的函数达到全局最优解，只有当函数为凸函数时，才能收敛到全局最优解。
2. 因为要在全数据集中优化参数，可能会导致时间太长。在海量训练数据下，这样做是十分耗时的。

梯度： 函数在该点处沿着该方向（此梯度的方向）变化最快，变化率最大。
对应于一元函数就是其一阶导数
对应于多元函数就是其一阶偏导数

梯度算法的计算公式为：$$w:=w+\alpha {▽}_{w}f(w)$$
其中，$\alpha$是步长，${▽}_{w}f(w)$是关于$w$的梯度。

2.二项Logistic回归模型参数估计

二项逻辑斯谛回归模型
$$P(Y=1|x)=\frac{exp(w\bullet x)}{1+exp(w\bullet x)}$$
$$P(Y=0|x)=\frac{1}{1+exp(w\bullet x)}$$
其中，$w$为权值向量，$x$为输入。
极大似然估计
设：$$\begin{gathered} P(Y=1|x)=\pi (x), \\ P(Y=0|x)=1-\pi (x) \\ {h}_w(x)=\frac{1}{1+exp(w\bullet x)} \end{gathered}$$
似然函数为：$$\prod _{i=1}^N[\pi (x^i{)}^{y^i}[1-\pi (x^i){]}^{1-y^i}$$
两边同时取对数，
$$\begin{array}{rl}L(w)& =\sum _{i=1}^N[y^{i}log\pi (x^i)+(1-y^i)log(1-\pi (x^i))]\\ & =\sum _{i=1}^N[y^{i}log\frac{\pi (x^i)}{1-\pi (x^i)}+log(1-\pi (x^i))]\\ & =\sum _{i=1}^N[y^i(w\bullet x^i)-log(1+exp(w\bullet x^i))]\end{array}$$
然后对每个权值向量$w_i求偏导，$
$$\begin{array}{rl}\frac{\partial L(w)}{\partial w_j}& =\sum _{i=1}^N{y}^i{x}_j^i-\sum _{i=1}^N\frac{exp(w\bullet x^i)}{1+exp(w\bullet x^i)}{x}_j^i\\ & =\sum _{i=1}^N[y^i-\frac{exp(w\bullet x^i)}{1+exp(w\bullet x^i)}]x_j^i\\ & =\sum _{i=1}^N[y^i-h_w(x^i)]x_j^i\end{array}$$
其中，$x_j^i$表示第$i$个样本的第$j$个特征。
由此，逻辑回归的最大似然估计的迭代公式为：$$w:=w+\sum _{i=1}^N[y^i-h_w(x^i)]x_j^i$$
注意：
关于Logistic回归的参数估计有其他推导版本，我在看的时候，在一个复合函数求导上纠结了半天，完整的展示如下：
$$h_{\theta }(x)=g({\theta }^{\tau }x)=\frac{1}{1+exp(-{\theta }^{\tau }x)}$$
一个样本的概率分布为：
$$P(Y=1|x;\theta )=h_{\theta }(x)P(Y=0|x;\theta )=1-h_{\theta }(x)$$
对回归系数$\theta$进行最大似然估计，
$$\prod _{i=1}^m{h}_{\theta }(x^i{)}^{y^i}(1-h_{\theta }(x^i){)}^{1-y^i}$$
取对数，
$$L(\theta )=\sum _{i=1}^m[y^{i}logh(x^i)+(1-y^i)log(1-h(x^i))]$$
求偏导(此处对复合函数求偏导，要小心)
$$\begin{array}{rl}\frac{\partial L(\theta )}{\partial {\theta }_j}& =\frac{y}{h({\theta }^{\tau }x)}\frac{\partial h({\theta }^{\tau }x)}{\partial {\theta }_j}-\frac{1-y}{1-h({\theta }^{\tau }x)}\frac{\partial h({\theta }^{\tau }x)}{\partial {\theta }_j}\\ & =[\frac{y}{h({\theta }^{\tau }x)}-\frac{1-y}{1-h({\theta }^{\tau }x)}]\frac{\partial h({\theta }^{\tau }x)}{\partial {\theta }_j}\\ & =[\frac{y}{h({\theta }^{\tau }x)}-\frac{1-y}{1-h({\theta }^{\tau }x)}]h({\theta }^{\tau }x)(1-h({\theta }^{\tau }x))\frac{\partial ({\theta }^{\tau }x)}{\partial {\theta }_j}\\ & =[y(1-h({\theta }^{\tau }x))-(1-y)h({\theta }^{\tau }x)]x_j\\ & =[y-h({\theta }^{\tau }x)]x_j\\ & =\sum _{i=1}^m[y^i-h({\theta }^{\tau }{x}^i)]x_j^i\end{array}$$
重点：第三行的 $h({\theta }^{\tau }x)(1-h({\theta }^{\tau }x)$ 怎么来的？带分数的指数函数求导，
$$h(z)=\frac{1}{1+e^{-z}}$$
等式两边同时取对数，
$$\begin{array}{rl}lnh(z)& =ln\frac{1}{1+e^{-z}}\\ & =-ln(1+e^{-z})\end{array}$$
等式两边同时关于$z$求导，
$$\frac{1}{h(z)}\frac{\partial h(z)}{\partial z}=-\frac{-e^{-z}}{(1+e^{-z})}$$
$$\frac{\partial h(z)}{\partial z}=\frac{e^{-z}}{(1+e^{-z})}h(z)$$
$$\begin{array}{rl}\frac{\partial h(z)}{\partial z}& =\frac{e^{-z}}{(1+e^{-z})}h(z)\\ & =\frac{(1+e^{-z})-1}{(1+e^{-z})}h(z)\\ & =h(z)(1-\frac{1}{(1+e^{-z})})\\ & =h(z)(1-h(z))\end{array}$$

3.随机梯度上升法

为了克服梯度上升法迭代速度慢的问题，提出来了随机梯度上升法。这个算法优化的不是在全部训练数据集上损失函数，而是在每一轮迭代中，随机优化某一条训练数据上的损失函数，这样每一轮参数的更新速度将大大加快。
之后在实战分析中，我们提出来改进的随机梯度上升算法，改进在

1. 学习率随着迭代次数会不断减小，但不会减小到0；
2. 通过随机选取样本来更新回归参数，这样会减少周期性波动；
3. 增加迭代次数

三、实战分析

1.梯度上升法

import numpy as np
import matplotlib.pyplot as plt
# 预处理数据
def loadDataSet():
    # 创建空列表，一个储存特征值，一个储存类别值
    dataMat = []; labelMat = []
    # 打开测试集
    fr = open('testSet.txt')
    # 读取测试集的每一行
    for line in fr.readlines():
        # 去掉每一行首尾的空格，并按照空格进行划分
        lineArr = line.strip().split()
        # 将测试集每一行的两个特征x1,x2,加上x0添加到dataMat中
        # 对x0的解释：原始形式 y=wx+b,为了统一表示，写成y=wx,此时，b作为了权重系数，其对应的特征值为1
        dataMat.append([1.0,float(lineArr[0]),float(lineArr[1])])
        # 将当前行的标签添加到labelMat中
        labelMat.append(int(lineArr[2]))
    return dataMat, labelMat

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

# 梯度上升算法
def gradAscent(dataMatIn,classLabels):
    # 将数据集列表转化成numpy矩阵
    dataMatrix = np.mat(dataMatIn)
    # 将标签列表转化成numpy矩阵，并进行转置
    labelMat = np.mat(classLabels).transpose()
    # 记录数据集的列数，即特征的个数
    n = np.shape(dataMatrix)[1]
    # 初始化步长
    alpha = 0.001
    # 初始化回归系数的权重
    weights = np.ones((n,1))
    # 最大的迭代次数
    MaxCircles = 500
    for i in range(MaxCircles):
        # 当前权重系数下的预测类别
        h = sigmoid(dataMatrix*weights)
        # 真实类别与预测类别的差
        error = labelMat-h
        # 权重系数的更新，后面的梯度被替换了，具体的推导公式见上文
        weights = weights + alpha*dataMatrix.transpose()*error
    return weights

判定边界： 在判定边界时，为什么会有$z=w_0{x}_0+w_1{x}_1+\cdots +w_n{x}_n=0$
我们观察Sigmoid函数，就会发现：
当$z>=0$时，$\sigma (z)>=0.5$,此时$y=1;$
当$z<=0$时，$\sigma (z)>=0.5$,此时$y=1;$
所以，$z=0$是一个决策边界，当它大于0或小于0时，逻辑回归模型分别预测不同的分类结果。

# 分析数据，画出决策边界
# 关于决策边界，上文有讲到
def plotBestFit(wei):
    # .getA()将括号里的数据集转化成数组
    weights = wei.getA()
    # 得到特征集和类别集
    dataMat,labelMat = loadDataSet()
    # 将dataMat转化成数组
    dataArr = np.array(dataMat)
    # 记录dataArr的行数
    n = np.shape(dataArr)[0]
    # 创建两个列表，储存类别1的后两个特征
    xcord1 = [];ycord1 = []
    # 创建两个列表，储存类别0的后两个特征
    xcord2 = [];ycord2 = []
    # 遍历数据集的每一行
    for i in range(n):
        # 如果为类别1
        if labelMat[i] == 1:
            # 将每一行得到的特征分别放在xcord1，ycord1中
            xcord1.append(dataArr[i,1]);ycord1.append(dataArr[i,2])
        else:
            # 否则，将类别0的一行得到的特征分别放在xcord1，ycord1中
            xcord2.append(dataArr[i,1]);ycord2.append(dataArr[i,2])
    # 创建画布
    fig = plt.figure()
    # 绘图区域，子图的个数
    ax = fig.add_subplot(111)
    # 散点图 s表示 marker 的大小；c表示 color;
    ax.scatter(xcord1,ycord1,s = 30,c = 'red',marker = 's')
    ax.scatter(xcord2,ycord2,s = 30,c = 'green',marker = 's')
    x=np.arange(-3.0,3.0,0.1)
    # 形如，w0+w1x1+w2x2=0,为什么是这个形式，见决策边界分析
    y=(-weights[0]-weights[1]*x)/weights[2]
    # 画直线
    plt.plot(x,y)
    # 加标签
    plt.xlabel('x1');plt.ylabel('x2')
    plt.show()

dataMat, labelMat = loadDataSet()# 将numpy矩阵转化成数组
weights = gradAscent(dataMat,labelMat)
print(plotBestFit(weights))

结果

2.随机梯度上升法

# 随机梯度上升
# 伪代码：
# 所有回归系数初始化为1
# 对数据集中的每个个体
#         计算该样本的梯度
#         使用alpha*gradient更新回归系数
# 返回回归系数值
def stocGradAscent0(dataMat,classlabels):
    dataMatrix = np.array(dataMat)
    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
        weights = weights + alpha*error*dataMatrix[i]
    return weights

weights = stocGradAscent0(dataMat,labelMat)
print(plotBestFit(weights))

注意：

def plotBestFit(wei):
    # .getA()将括号里的数据集转化成数组
    # 这一行在梯度上升中将矩阵转化成数组，这里不用转化了
    #weights = wei.getA()
    weights = wei
    # 得到特征集和类别集
    dataMat,labelMat = loadDataSet()
    # 将dataMat转化成数组
    dataArr = np.array(dataMat)

结果：

效果可能没有梯度上升法好，但这个结果是很正常的，因为梯度上升法是在整个数据集上迭代了500次才得到的。

3.改进的随机梯度上升法

# 改进的随机梯度上升法
def stocGradAscent1(dataMat,classlabels):
    # 将数据集转化成numpy数组
    dataMat = np.array(dataMat)
    # 得到数据集的行和列
    m,n = np.shape(dataMat)
    # 最大的迭代次数
    numIter = 200
    # 初始化权重向量
    weights = np.ones(n)
    # 开始进行迭代
    for i in range(numIter):
        # 获取数据集的下标列表
        dataIndex = list(range(m))
        # 遍历行列表
        for j in range(m):
            # 更新步长
            alpha = 4/(1.0+i+j)+0.01
            # 随机获取样本对权重进行更新
            randomIndex = int(np.random.uniform(0, len(dataIndex)))
            # 计算当前实例对应的Sigmoid值
            h = sigmoid(sum(dataMat[randomIndex]*weights))
            # 计算当前样本的残差
            error = classlabels[randomIndex] - h
            # 更新权值，注意和之前不一样的地方是：用的是向量的乘法，即对应位置相乘
            weights = weights + alpha * error *dataMat[randomIndex]
            # 删除数据集中已经被选的实例下标
            del (dataIndex[randomIndex])
    return weights

原因：python 2返回的是列表，python 3返回的是迭代器，range 返回的是range 对象
办法: 将range() 改为 list(range())
结果：

weights = stocGradAscent1(dataMat,labelMat)
print(plotBestFit(weights))# 输出分类图

4.实例：从疝气病症预测病马的死亡率

预测过程：

1. 收集数据：给定数据文件
2. 准备数据：用Python解析文本文件并填充缺失值
3. 分析数据：可视化并观察数据
4. 训练算法：使用优化算法，确定最佳的系数
5. 测试算法：为了量化回归的效果，需要观察错误率，根据错误率决定是否退回到训练阶段，通过改变迭代次数和步长等参数来得到更好的回归系数
6. 使用算法：输入马的症状，得到预测结果

准备数据：处理数据中的缺失值
数据中的缺失值如何处理一直是一个重要的问题，比如，假设有100个样本和20个特征，这些数据都是机器收集的，但如果因为机器上的某一个传感器失效，导致某一个特征失效怎么办呢？放弃整个数据集显然是不可取的。
一般而言，有如下的处理方法：
1） 当缺失值较多的特征处理
一般如果某特征的缺失量过大，我们会直接将该特征舍弃掉，否则可能反倒会带入较大的noise，对结果造成不良影响。
2）缺失值较少的特征处理（重点）

1. 把NaN直接作为一个特征，假设用0表示
2. 用特征的均值来填充缺失值，将所有行用各自的均值填充
3. 也可以指定某些行进行填充 ，训练集中有缺失值，而测试集中无缺失值。这时候最适合的处理方式应该是对缺失值取条件均值或条件中值，即即根据该用户的label值类别，取所有该label下用户该属性的均值或中值。
4. 用上下数据进行填充；
5. 用插值法填充，.interpolate()；
6. 预测模型（Prediction Model）,预测模型是处理缺失值的复杂方法之一， 通过创建一个预测模型来估计替代缺失值。 在这种情况下，我们将数据集分为两组：一组没有变量的缺失值，另一组有缺少值， 第一个数据集成为模型的训练数据集，而具有缺失值的第二个数据集是测试数据集，变量与缺失值被视为目标变量。 接下来，我们创建一个模型，根据训练数据集的其他属性预测目标变量，并填充测试数据集的缺失值；
7. 使用KNN算法等预测缺失值

算法代码：

# 分类函数
def classifyVector(inX,weights):
    # 计算此时实例对应的Sigmoid值
    prob = sigmoid(sum(inX*weights))
    # 如果大于 0.5，则属于类1；如果小于0.5，则对应类0
    if prob > 0.5 : return 1.0
    else :   return 0.0
def colicTest():
    # 打开训练数据集
    frTrain = open('horseColicTraining.txt')
    # 打开测试数据集
    frTest = open('horseColicTest.txt')
    # 创建两个空列表来储存训练数据集和标签集
    trainingSet = [];trainingLabel = []
    # 遍历数据集的每一行
    for line in frTrain.readlines():
        # 对当前行按制表符进行分割
        currline = line.strip().split('\t')
        # 新建列表用来储存每个样本的特征向量
        lineArr = []
        # 将每一个样本的0-20列储存在lineArr中，即特征列表
        for  i in range(21):
            lineArr.append(float(currline[i]))
        trainingSet.append(lineArr)
        # 将每一个样本的最后一列储存在trainLabel中
        trainingLabel.append(float(currline[21]))
    # 调用随机梯度上升法更新logistic回归的权值参数
    trainingWeights = stocGradAscent1(trainingSet,trainingLabel)
    # 初始化误分类的个数和训练集的个数为0
    errorCount = 0; numTestVec = 0
    # 遍历测试集的每一行
    for line in frTest.readlines():
        # 记录测试集个数
        numTestVec +=1.0
        # 对当前样本按照制表符进行分割
        currline = line.strip().split('\t')
        # 新建列表用来储存每个样本的特征向量
        lineArr = []
        # 将每一个样本的0-20列储存在lineArr中，即特征列表
        for i in range(21):
            lineArr.append(float(currline[i]))
        # 利用分类预测函数对该样本进行预测，并与样本标签进行比较
        if int(classifyVector(lineArr,trainingWeights)) !=int(currline[21]):
            # 如果误分类，则计数器加1
            errorCount +=1
    # 计算误分类率
    errorate =(float(errorCount)/numTestVec)
    print('The error rate of the test is %f:' % errorate)
    return errorate
# 计算10次分类结果的均值
def multiTest():
    # 迭代10次，初始化误分类的和为0
    numTests=10;errorSum = 0.0
    # 开始迭代
    for i in range(numTests):
        errorSum +=colicTest()
    # 10次误分类的总数/10
    ave_error = errorSum/float(numTests)
    print('After %d iterations the average error rate is :%f'%(numTests,ave_error))
    print('After %d iterations the average right rate is :%f'%(numTests,1-ave_error))

结果：

print(multiTest())

资料参考：
1.《统计学习方法》李航
2.《机器学习实战》Peter Harrington
3.逻辑回归（logistic regression）原理详解：https://blog.csdn.net/guoziqing506/article/details/81328402
4.机器学习之梯度与梯度下降法：https://blog.csdn.net/qq_20412595/article/details/81409744
5.梯度下降法（上升法）的几何解释：https://blog.csdn.net/xmu_jupiter/article/details/22220987
6.常用数学符号的 LaTeX 表示方法:http://mohu.org/info/symbols/symbols.htm
7.逻辑回归（Logistic Regression, LR）简介:https://blog.csdn.net/jk123vip/article/details/80591619
8.机器学习系列(1)_逻辑回归初步：https://blog.csdn.net/longxinchen_ml/article/details/49130931