Logistic Regression

《Logistic Regression》

本文源代码以及数据均可从:点击此处获取

1. Logistic Regression 模型：

1.1 线性可分 VS 线性不可分：

对于一个分类问题，通常可以分为线性可分与线性不可分两种。如果一个分类问题可以使用线性判别函数正确分类，则称该问题为线性可分；否则为线性不可分问题。

图1.1 线性可分

图1.2 线性不可分

图1.3 线性不可分转化为线性可分

1.2 Logistic Regression 模型

Logistic Regression 模型是广义线性模型的一种，属于线性的分类模型。对于图1.1所示的线性可分的问题，需要找到一条直线，能够将两个不同的类区分开，这条红色的直线被称为**超平面**。

对于上述的超平面，可以使用如下的线性函数表示：

$$Wx+b=0$$

其中W为权重，b为偏置。若在多为的情况下，权重W和偏置b均为向量。在Logistic Regression算法中，通过对训练样本的学习，最终得到该超平面，将数据分成正负两个类别。此时，可以使用阈值函数，将样本映射到不同的类别中，常见的阈值函数有Sigmoid函数，其形式如下所示：

$$f(x)=\frac{1}{1+e^{-x}}$$

Sigmoid函数的图像如下图所示：

图1.2 Sigmoid 函数

从Sigmoid函数的图像上可以看出，其函数的值域为（0，1），在0附近的变化比较明显。其导函数$f^{\prime }(x)$为：

$$f^{\prime }(x)=\frac{e^{-x}}{(1+e^{-x}{)}^2}=f(x)[1-f(x)]$$

现在，让我们利用Python实现Sigmoid函数，为了能够使用numpy中的函数，我们首先需要导入 numpy :

import numpy as np

def sig(x):
    """Sigmoid函数
    Input: x(mat) : feature * w
    Output: sigmoid(x)(mat):Sigmoid值
    """
    return 1.0 / (1 + np.exp(-x))

上述程序中，Sigmoid函数的输出为Simoid值。对于输入向量X，其属于正例的概率为：

$$P(y=1|X,W,b)=\sigma (WX+b)=\frac{1}{1+e^{-(WX+b)}}$$

其中，$\sigma$表示的是Sigmoid函数。那么，对于输入向量 X ，其属于负例的概率为：

$$P(y=1|X,W,b)=1-P(y=1|X,W,b)=1-\sigma (WX+b)=\frac{e^{-(WX+b)}}{1+e^{-(WX+b)}}$$

对于Logistic Regression 算法来说，需要求解的分隔超平面中的参数，即为权重矩阵 W 和偏置向量 b ,那么，这些参数该如何求解呢？为了求解模型的两个参数，首先必须定义损失函数。

1.3 损失函数：

对于上述的Logistic Regression算法，其属于类别y的概率为：

$$P(y|X,W,b)=\sigma (WX+b{)}^y(1-\sigma (WX+b){)}^{1-y}$$

要求上述问题中的参数W和b，可以使用极大似然法对其进行估计。假设训练数据集有m个训练样本$\left\{(X^{(1)},y^{(1)}),(X^{(2)},y^{(2)}),\cdots ,(X^{(m)},y^{(m)})\right\}$,则其似然函数为：

$$L_{W,b}=\prod _{i=1}^m\left[h_{W,b}(X^{(i)}{)}^{y^{(i)}}(1-h_{W,b}(X^{(i)}){)}^{1-y^{(i)}}\right]$$

其中，假设函数$h_{W,b}(X^{(i)})$为：

$$h_{W,b}(X^{(i)})=\sigma (W{X}^{(i)}+b)$$

对于似然函数的极大值的求解，通常使用Log似然函数，在Logistic Regression算法中，通常是将负的Log似然函数作为其损失函数，即 the negative log-likelihood(NLL)作为其损失函数，此时，需要计算的是NLL的极小值。损失函数$l_{W,b}$为：

$$l_{W,b}=-\frac{1}{m}log{L}_{W,b}=-\frac{1}{m}\sum _{i=1}^m\left[y^{(i)}log(h_{W,b}(X^{(i)}))+(1-y^{(i)})log(1-h_{W,b}(X^{(i)}))\right]$$

此时，我们需要求解的问题为：

$$\underset{W,b}{\min }{l}_{W,b}$$

为了求得损失函数$l_{W,b}$的最小值，可以使用基于梯度的方法进行求解。

2. 梯度下降法：

在机器学习算法中，对于很多监督学习模型，需要对原始的模型构建损失函数l，接下来表示通过优化算法对损失函数l进行优化，以便寻找到最优的参数W。在求解机器学习参数W的优化算法时，使用较多的是基于梯度下降的优化算法（Gradient Descent,GD）。

梯度下降法有很多优点，其中，在梯度下降法的求解过程中，只需求解损失函数的一阶导数，计算的成本比较小，这使得梯度下降法能在很多大规模数据集上得到应用。梯度下降法的含义是通过当前点的梯度方向寻找到新的迭代点，并从当前点移动到新的迭代点继续寻找新的迭代点，直到找到最优解。

2.1 梯度下降法的流程：

梯度下降法是一种迭代型的优化算法，根据初始点在每一次迭代的过程中选择下降方向，进而改变需要修改的参数，对于优化问题min f(w),梯度下降法的详细过程如下所示：

• 随机选择一个初始点$w_0$
• 重复以下过程：
  ° 决定梯度下降的方向： $d_i=-\frac{\partial }{\partial w}f(w){|}_{w_i}$
  ° 选择步长: $\alpha$
  ° 更新: $w_{i+1}=w_i+\alpha d_i$
• 直到满足终止条件

图1.4 梯度下降的过程

在初始时，在点$w_0$处，选择下降的方向$d_0$，选择步长 $\alpha$，更新w的值，此时到达$w_1$处，判断是否满足终止的条件，发现并未到达最优解$w^{\ast }$，重复上述的过程，直至到达$w^{\ast }$。

2.2 凸优化与非凸优化：

简单来说，凸优化问题是指只存在一个最优解的优化问题，即任何一个局部最优解即是全局最优解，如图1.5所示。

图1.5 凸函数

非凸优化是指解空间中存在多个局部最优解，而全局最优解是其中的某一个局部最优解，如图1.6所示。

图1.6 非凸函数

最小二乘（Least Squares）、岭回归（Ridge Regression）和逻辑回归（Logistic Regression）的损失函数都是凸优化问题。

2.3 利用梯度下降法训练Logistic Regression模型：

对于上述的Logistic Regression算法的损失函数可以通过梯度下降法对其进行求解，其梯度为：

$${▽}_{W_j}(l_{W,b})=-\frac{1}{m}\sum _{i=1}^m(y^{(i)}-h_{W,b}(X^{(i)}))x_j^{(i)}$$

$${▽}_b(l_{W,b})=-\frac{1}{m}\sum _{i=1}^m(y^{(i)}-h_{W,b}(X^{(i)}))$$

其中，$x_j^{(i)}$表示的是样本$X^{(i)}$的第j个分量。取$w_0=b$,且将偏置项的变量$x_0$设置为1，则可以将上述的梯度合并为：

$${▽}_{W_j}(l_{W,b})=-\frac{1}{m}\sum _{i=1}^m(y^{(i)}-h_{W,b}(X^{(i)}))x_j^{(i)}$$

根据梯度下降法，得到如下的更新公式：

$$W_j=W_j+\alpha {▽}_{W_j}(l_{W,b})$$

利用上述的Logistic Regression中权重的更新公式，我们可以实现Logistic Regression模型的训练，利用梯度下降法训练模型的具体过程，如下程序所示：

# 程序清单 1-2 Logistic Regression 模型的训练

def lr_train_bgd(feature, label, maxCycle, alpha):
    """利用梯度下降法训练LR模型
    Input: feature(mat)特征
           label(mat) 标签
           maxCycle(int)最大迭代次数
           alpha(float)学习率
    Output: w(mat)权重
    """
    n = np.shape(feature)[1]                   #特征个数
    w = np.mat(np.ones((n, 1)))                #初始化权重
    i = 0
    while i <= maxCycle:                       # 在最大迭代次数的范围内
        i += 1                                 #当前的迭代次数
        h = sig(feature * w)                   # 计算Sigmoid值
        err = label - h
        if i % 100 == 0:
            print("\t----------iter=" + str(i) + " , train error rate= " + str(error_rate(h, label)))
        w = w + alpha * feature.T * err    # 权重修正
    return w
    

在上述程序清单中，函数lr_train_bgd使用了梯度下降法对Logistic Regression算法中的损失函数进行优化，其中，lr_train_bgd函数的输入为训练样本的特征、训练样本的标签、最大的迭代次数和学习率，在每一次迭代的过程中，需要计算当前的模型的误差，误差函数为 error_rate，error_rate函数的具体实现形式如程序清单1-3所示。在迭代的过程中，不断通过梯度下降的方法对Logistic Regression算法中的权重进行更新。

#程序清单 1-3 error_rate 函数的实现

def error_rate(h, label):
    '''计算当前的损失函数值
    Input: h(mat):预测值
           label(mat): 实际值
    Output: err/m(float): 错误率
    '''
    
    m = np.shape(h)[0]
    
    sum_err = 0.0
    for i in range(m):
        if h[i, 0] > 0 and (1 - h[i, 0]) > 0:
            sum_err -= (label[i,0] * np.log(h[i,0]) + (1 - label[i,0])*np.log(1 - h[i,0] ))
        else:
            sum_err -= 0
    return sum_err / m

3. 梯度下降法的若干问题：

3.1 选择下降的方向：

为了求解优化问题f(w)的最小值，我们希望每次迭代的结果能够接近最优值 $w^{\ast }$,对于一维的情况，如图 1.7 所示：

图1.7 下降的方向

若当前点的梯度为负，则最小值在当前点的右侧，若当前点的梯度为正，则最小值在当前点的左侧，负的梯度即为下降的方向。对于上述的一维的情况，有下述的更新规则：

$$w_{i+1}=w_i-{\alpha }_i\frac{df}{dw}{|}_{w=w_i}$$

其中， ${\alpha }_i$为步长。对于二维的情况，此时更新的规则如下：

$$W_{i+1}=W_i-{\alpha }_i▽f(W_i)$$

3.2 步长的选择：

对于步长 $\alpha$ 的选择，若选择太小，会导致收敛的速度比较慢；若选择太大，则会出现震荡的现象，即跳过最优解，在最优解附近徘徊，上述两种情况如图1.8 所示：

图1.8 步长太长或者太小

因此，选择合适的步长对于梯度下降法的收敛效果显得尤为重要，如图 1.9 所示：

图1.9 合适的步长

4. Logistic Regression 算法实践：

有了以上的理论准备，接下来，我们利用已经完成的函数，构建Logistic Regression 分类器。我们利用线性可分的数据作为训练样本训练Logistic Regression 模型，在构建模型的过程中，主要分为两个步骤：
（1）利用训练样本训练模型；
（2）利用训练好的模型对新样本进行预测。

4.1 利用训练样本训练 Logistic Regression 模型：

首先，我们利用训练样本训练模型，为了使得Python能够支持中文的注释和利用numpy工具，我们需要在训练文件的开始加入：

# coding:UTF-8
import numpy as np

在训练模型中，其主函数如程序清单1-4 所示：

# 程序清单 1-4 训练模型的主函数

if __name__=="__main__":
    # 1. 导入训练数据
    print("-------- 1. load data --------")
    feature, label = load_data1("data.txt")
    
    # 2. 训练 LR 模型
    print("-------- 2. training --------")
    w = lr_train_bgd(feature, label, 1000, 0.01)
    
    # 3. 保存最终的模型
    print("-------- 3. save model --------")
    save_model("weights", w)
    

-------- 1. load data --------
-------- 2. training --------
	----------iter=100 , train error rate= 0.0011343552118725198
	----------iter=200 , train error rate= 0.0009477843077847437
	----------iter=300 , train error rate= 0.0008150655965653146
	----------iter=400 , train error rate= 0.0007156807636573233
	----------iter=500 , train error rate= 0.0006383910251337722
	----------iter=600 , train error rate= 0.0005765152961049185
	----------iter=700 , train error rate= 0.0005258283298582945
	----------iter=800 , train error rate= 0.00048352514286009306
	----------iter=900 , train error rate= 0.0004476693385107398
	----------iter=1000 , train error rate= 0.0004168803688943547
-------- 3. save model --------

#  程序清单 1-5 导入训练数据的 load_data 函数
import numpy as np

def load_data1(file_name):
    '''
    Input: file_name(string) 训练数据的位置
    Output: feature_data(mat)特征
            label_data(mat)标签
    '''
    f = open(file_name) #打开文件
    feature_data = []
    label_data = []
    for line in f.readlines():
        feature_tmp = []
        label_tmp = []
        lines = line.strip().split("\t")
        feature_tmp.append(1)  # 偏置项
        for i in range(len(lines) - 1):
            feature_tmp.append(float(lines[i]))
        label_tmp.append(float(lines[-1]))
        
        feature_data.append(feature_tmp)
        label_data.append(label_tmp)
    f.close()  # 关闭文件
    return np.mat(feature_data), np.mat(label_data)

在 load_data 函数中，其输入的为训练数据所在的位置，其输出的为训练数据的特征和训练数据的标签。

# 程序清单 1-6 保存最终的模型的 save_model 函数
def save_model(file_name, w):
    '''
    保存最终的模型
    Input: file_name(string):模型保存的文件名
           w(mat): LR 模型的权重
    '''
    
    m = np.shape(w)[0]
    f_w = open(file_name, "w")
    w_array = []
    for i in range(m):
        w_array.append(str(w[i, 0]))
    f_w.write("\t".join(w_array))
    f_w.close()

在程序清单1-6中，函数save_model将训练好的LR模型以文件的形式保存，其中 save_model函数的输入为保存的文件名file_name和所需保存的模型w。

数据可视化：

import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
feature, label = load_data("data.txt")
w = np.mat(np.ones((3, 1)))

#h = sig(feature * w)
#print(feature)
#n = np.shape(feature)[1]
#print(n)

h = sig(feature * w)
print(np.shape(h))

(200, 1)

4.2 对新数据进行预测：

对于分类算法而言，训练好的模型需要能够对新的数据集进行划分。利用上述步骤，我们训练好LR模型，并将其保存在“weights”文件中，此时，我们需要利用训练好的LR模型对新数据进行预测，同样，为了能够使用numpy中的函数和对中文注释的支持，在文件“lr_test.py”开始，我们加入：

# coding: UTF-8
import numpy as np

主函数如下程序清单1-7所示。

# 程序清单1-7 测试的主函数

if __name__=="__main__":
    # 1. 导入 LR 模型
    print("---------- 1. load model ----------")
    w = load_weight("weights")
    n = np.shape(w)[1]
    
    # 2. 导入测试数据
    print("---------- 2. load test data ----------")
    testData = load_data("test_data", n)
    
    # 3. 对测试数据进行预测
    print("---------- 3. get prediction ----------")
    h = predict(testData, w)
    
    # 4. 保存最终预测结果
    print("---------- 4. save prediction ----------")
    save_result("result", h)

---------- 1. load model ----------
---------- 2. load test data ----------
---------- 3. get prediction ----------
---------- 4. save prediction ----------

# 程序清单 1-8 导入模型的 load_weight 函数
def load_weight(w):
    '''导入LR模型
    Input: w(string):权重所在的文件位置
    Output: np.mat(w)(mat)权重的矩阵
    '''
    f = open(w)
    w = []
    for line in f.readlines():
        lines = line.strip().split("\t")
        w_tmp = []
        for x in lines:
            w_tmp.append(float(x))
        w.append(w_tmp)
    f.close()
    return np.mat(w)

在程序清单1-8中，首先需要导入numpy模块和lr_train中的sig函数。在load_weight函数中，其输入是权重所在的位置，在导入函数中，将其数值导入到权重矩阵中。

#程序清单 1-9 导入测试集的 load_data 函数
def load_data(file_name, n):
    '''导入测试数据
    Input: file_name(string):测试集的位置
           n (int):特征的个数
    Output: np.mat(feature_data) (mat):测试集的特征
    '''
    f = open(file_name)
    feature_data = []
    for line in f.readlines():
        feature_tmp = []
        lines = line.strip().split("\t")
        # print(lines[2])
        if len(lines) != n-1:    #特征的个数用于判断测试集是否符合要求，若不符合要求，则丢弃。
            continue
        feature_tmp.append(1)
        for x in lines:
            # print(x)
            feature_tmp.append(float(x))
        feature_data.append(feature_tmp)
    f.close()
    return np.mat(feature_data)

在导入测试集的 load_data 函数中，其输入为测试集的位置和特征的个数，其中特征的个数用于判断测试集是否符合要求，若不符合要求，则丢弃。

# 程序清单 1-10 对新数据集进行预测的 predict 函数
def predict(data, w):
    '''对测试数据进行预测
    Input: data(mat): 测试数据的特征
           w(mat): 模型的参数
    Output: h(mat): 最终的预测结果
    '''
    h = sig(data * w.T) # 取得 Sigmoid 值
    m = np.shape(h)[0]
    for i in range(m):
        if h[i, 0] < 0.5:
            h[i, 0] = 0.0
        else:
            h[i, 0] = 1.0
    return h

在 predict 函数中，其输入为测试数据的特征和模型的权重，输出为最终的预测结果。通过特征与权重的乘积，再对其求 Sigmoid 函数值得到最终的预测结果。在此，使用到了文件“lr_train.py”中的 sig 函数，因此，在文件“lr_test.py”中，需要导入 sig 函数：

from lr_train import sig

在计算最终的输出时，为了将Sigmoid函数输出的概率值转换成{0,1}，通常可以取0.5作为边界。

# 程序清单 1-11 保存最终预测结果的 save_result 函数
def save_result(file_name, result):
    '''保存最终的预测结果
    Input:  file_name(string):预测结果保存的文件名
            result(mat): 预测的结果
    '''
    m = np.shape(result)[0]
    #输出预测结果到文件
    tmp = []
    for i in range(m):
        tmp.append(str(result[i, 0]))
    f_result = open(file_name, "w")
    f_result.write("\t".join(tmp))
    f_result.close()

在程序清单 1-11 中，函数 save_result 实现将预测结果存到指定的文件中，函数 save_result 的输入为预测结果保存的文件名 file_name 和预测的结果 result，最终将 result 中的数据写入到文件 file_name 中。

参考文献

---

[1] 李航.《统计学习方法》[M].北京.清华大学出版社.2012.
[2] 周志华.《机器学习》[M].北京.清华大学出版社.2016.
[3] Peter Harrington. 机器学习实战[M].王斌，译.北京.人民邮电出版社.2013.
[4] Chapelle O, Manavoglu E, Rosales R. Simple and Scalable Response Prediction for Display Advertising[J]. Acm Transactions on Intelligent Systems & Technology,2014,5(4):1-34.