2017-08-23

逻辑回归

文章目录

原文地址:https://ngunauj.github.io

逻辑回归（Logistic Regression）是机器学习中的一种分类模型，由于算法的简单和高效，在实际中应用非常广泛。其实质实在线性回归的基础上，套用了一个逻辑函数，但也就由于这个逻辑函数，使得逻辑回归模型成为了机器学习领域一颗耀眼的明星，更是计算广告学的核心。通常逻辑回归用在预测一个用户是否点击特定的商品,判断用户的性别,预测用户是否会购买给定的品类,判断一条评论是正面的还是负面的…本文主要记录下逻辑回归的基础知识，更多的内容以后再去搞吧

逻辑回归模型
最简单的回归是线性回归，在此借用Andrew NG的讲义，有如图1.a所示，X为数据点——肿瘤的大小，Y为观测值——是否是恶性肿瘤。通过构建线性回归模型，如hθ(x)所示，构建线性回归模型后，即可以根据肿瘤大小，预测是否为恶性肿瘤hθ(x)≥.05为恶性，hθ(x)<0.5为良性。

然而线性回归的鲁棒性很差，例如在图1.b的数据集上建立回归，因最右边噪点的存在，使回归模型在训练集上表现都很差。这主要是由于线性回归在整个实数域内敏感度一致，而分类范围，需要在[0,1]。逻辑回归就是一种减小预测范围，将预测值限定为[0,1]间的一种回归模型，其回归方程与回归曲线如图2所示。逻辑曲线在z=0时，十分敏感，在z>>0或z<<0处，都不敏感，将预测值限定为(0,1)。

Sigmoid 函数
说到逻辑回归，Sigmoid是一大要点，其表达式为：
$$g(z) = \frac{1}{1 + e^{-z} } \tag{1} \label{1}$$
它是一个可导函数，定义域为(−∞,+∞)(−∞,+∞)，值域为[0, 1]，其导数为：
$$g’(z) = g(z)(1-g(z)) \tag{2} \label{2}$$说明一下，表达式(1)等价于使用线性回归模型的预测结果直接去逼近真实标记的对数几率，因此将其称作“对数几率回归（logit regression）”。使用这种方法有以下几大优点：

直接对样本进行建模，无需对样本进行先验假设；
其结果不仅可以预测出“label”，还可以得到近似的概率预测值；
sigmoid函数的数学性质良好，它是任意阶可导的凸函数，因此许多的优化方法都可使用。

极大似然估计(MLE)与损失函数
在机器学习理论中，损失函数（loss function）是用来估量你模型的预测值f(x)与真实值Y的不一致程度，它是一个非负实值函数，通常使用L(Y,f(x))来表示，损失函数越小，模型的鲁棒性就越好。损失函数是经验风险函数的核心部分，也是结构风险函数重要组成部分。模型的结构风险函数包括了经验风险项和正则项，通常可以表示成如下式子：
$$\theta^*=\arg \min_\theta \frac{1}{N}{}\sum_{i=1}^{N} L(y_i, f(x_i; \theta)) + \lambda\ \Phi(\theta) \tag{3}\label{3}$$
对于逻辑回归，其loss function是log损失，可以通过极大似然估计进行推导得到。首先，给定一个样本x，可以使用一个线性函数对自变量进行线性组合，
$$\theta_0 + \theta_1 x_1 + \theta_2 x_2 + \cdots + \theta_nx_n =\theta^Tx \tag{4} \label{4}$$
根据sigmoid函数，我们可以得出预测函数的表达式：$$h_{\theta}(x) = g(\theta^Tx) = \frac{1}{1 + e^{-\theta^Tx}} \tag{5} \label{5}$$
式(4)表示y=1时预测函数为hθ(x)。在这里，假设因变量y服从伯努利分布，取值为0和1，那么可以得到下列两个式子：
$$p(y=1 | x) = h_{\theta} (x) \tag{6} \label{6}$$
$$p(y=0 | x) = 1 - h_{\theta} (x) \tag{7} \label{7}$$
而对于上面的两个表达式，通过观察，我们发现，可以将其合并为表达式(8)(概率论中有讲)：
$$p(y | x) = h_{\theta} (x)^y (1-h_{\theta} (x))^{1-y} \tag{8} \label{8}$$
根据上面的式子，给定一定的样本之后，我们可以构造出似然函数，然后可以使用极大似然估计MLE的思想来求解参数。
但是，为了满足最小化风险理论，我们可以将MLE的思想转化为最小化风险化理论，最大化似然函数其实就等价于最小化负的似然函数。对于MLE，就是利用已知的样本分布，找到最有可能（即最大概率）导致这种分布的参数值；或者说是什么样的参数才能使我们观测到目前这组数据的概率最大。使用MLE推导LR的loss function的过程如下。
首先，根据上面的假设，写出相应的极大似然函数（假定有m个样本）：
$$\begin{aligned} L(\theta) &= \prod_{i=1}^{m} p(y^{(i)} | x^{(i)}; \theta) \\ &= \prod_{i=1}^{m} h_{\theta} (x^{(i)})^{y^{(i)}} (1-h_{\theta} (x^{(i)}))^{1-y^{(i)}} \ \end{aligned} \tag{9} \label{9}$$
直接对上面的式子求导会不方便，因此，为了便于计算，我们可以对似然函数取对数，经过化简可以得到下式的推导结果：
$$\begin{aligned} \log L(\theta) &= \sum_{i=1}^{m} \log \left [ (h_{\theta} (x_i)^{y^{(i)}} (1-h_{\theta} (x^{(i)}))^{1-y^{(i)}}) \right ] \\ &= \sum_{i=1}^{m} \left [ y^{(i)} \log h_{\theta} (x^{(i)}) + (1-y^{(i)}) \log(1-h_{\theta} (x^{(i)})) \right ] \ \end{aligned} \tag{10} \label{10}$$
因此，损失函数可以通过最小化负的似然函数得到，即下式：
$$J(\theta) = - \frac{1}{m} \sum_{i=1}^m \left [ y^{(i)} \log h_{\theta}(x^{(i)}) + (1-y^{(i)}) \log(1-h_{\theta}(x^{(i)})) \right ] \tag{11} \label{11}$$
在有的资料上，还有另一种损失函数的表达形式，但本质是一样的，如下：
$$J(\theta) = \frac{1}{m} \sum_{i=1}^m log(1 + e^{-y^{(i)} \theta^T x})$$
…之所以有人认为逻辑回归是平方损失，是因为在使用梯度下降来求最优解的时候，它的迭代式子与平方损失求导后的式子非常相似，从而给人一种直观上的错觉。

Gradient descent
梯度下降法又叫做最速下降法，为了求解使损失函数J(θ)最小时的参数θ，这里就以梯度下降为例进行求解，其迭代公式的推导过程如下：
$$\begin{aligned} \frac{ \partial J(\theta)} {\partial \theta_j} &= -\frac{1}{m} \sum_{i}^{m} \left [ y^{(i)}(1 - h_{\theta}(x^{(i)})) \cdot (-x_{j}^{(i)}) + (1 - y^{(i)}) h_{\theta} (x^{(i)}) \cdot (x_{j}^{(i)}) \right ] \\ &= - \frac{1}{m} \sum_{i}^{m} (-y^{(i)} \cdot x_j^{(i)} + h_{\theta}(x^{(i)}) \cdot x_j^{(i)}) \\ &= -\frac{1}{m} \sum_{i}^{m} (h_{\theta}(x^{(i)}) - y^{(i)}) x_j^{(i)} \end{aligned} \tag{12} \label{12}$$
通过上面得到，可以得到最后的迭代式子：
$$\theta_j = \theta_j - \alpha \sum_{i=1}^{m} (h_{\theta}(x^{(i)}） - y^{(i)}) x_j^{(i)} \tag{13}\label{13}$$
其中α是步长。求导过程 reference
最优化算法并不限于梯度下降，还有：

Newton Method（牛顿法）
Conjugate gradient method(共轭梯度法)
Quasi-Newton Method(拟牛顿法)
BFGS Method
L-BFGS(Limited-memory BFGS)

Regulization
上面提到了，结构风险函数包括了经验风险项和正则项，加入正则项相当于对参数加入了一个先验分布，常用的有L1和L2正则，L1，L2正则化项对模型的参数向量进行“惩罚”，从而避免单纯最小二乘问题的过拟合问题。正则化项本质上是一种先验信息，整个最优化问题从贝叶斯观点来看是一种贝叶斯最大后验估计，其中正则化项对应后验估计中的先验信息，损失函数对应后验估计中的似然函数，两者的乘积即对应贝叶斯最大后验估计的形式，如果将这个贝叶斯最大后验估计的形式取对数，即进行极大似然估计，你就会发现问题立马变成了损失函数+正则化项的最优化问题形式。
在逻辑回归求解中，如果对样本加上一个先验的服从高斯分布的假设，那么就取log后，式子经过化简后就变成在经验风险项后面加上一个正则项，此时损失函数变为。
$$J(\theta) = - \frac{1}{m}\sum_{i=1}^m \left [ y^{(i)} \log h_{\theta}(x^{(i)}) + (1-y^{(i)}) \log(1-h_{\theta}(x^{(i)})) \right ] + \frac{\lambda}{2m} \sum_{j=1}^{n} \theta_{j}^2 \tag{14}\label{14}$$
关于LR，这里还有个PDF可以参考一下：Lecture 6: logistic regression.pdf

模型·评估
对于LR分类模型的评估，常用AUC来评估，关于AUC的更多定义与介绍，可见参考文献2，在此只介绍一种极简单的计算与理解方法。对于下图的分类：

对于训练集的分类，训练方法1和训练方法2分类正确率都为80%，但明显可以感觉到训练方法1要比训练方法2好。因为训练方法1中，5和6两数据分类错误，但这两个数据位于分类面附近，而训练方法2中，将10和1两个数据分类错误，但这两个数据均离分类面较远。
AUC正是衡量分类正确度的方法，将训练集中的label看两类{0，1}的分类问题，分类目标是将预测结果尽量将两者分开。将每个0和1看成一个pair关系，团中的训练集共有5*5=25个pair关系，只有将所有pair关系一至时，分类结果才是最好的，而auc为1。在训练方法1中，与10相关的pair关系完全正确，同样9、8、7的pair关系也完全正确，但对于6，其pair关系(6，5)关系错误，而与4、3、2、1的关系正确，故其auc为(25-1)/25=0.96；对于分类方法2，其6、7、8、9的pair关系，均有一个错误，即(6,1)、(7,1)、(8,1)、(9,1)，对于数据点10，其正任何数据点的pair关系，都错误，即(10,1)、(10,2)、(10,3)、(10,4)、(10,5)，故方法2的auc为(25-4-5)/25=0.64，因而正如直观所见，分类方法1要优于分类方法2。

代码实现
我们使用Iris数据集数据集数据中特征向量有4个，分3个类别(数据包含三种鸢尾花的四个特征，分别是花萼长度(cm)、花萼宽度(cm)、花瓣长度(cm)、花瓣宽度(cm))。但是，我们为了数据的可视化，我们只保留2个特征花瓣长度(petal),花萼长度(sepal)

# -*- coding: utf-8 -*-
"""
Spyder Editor
This is a temporary script file.
"""
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np
df = pd.read_csv('http://archive.ics.uci.edu/ml/machine-learning-databases/iris/iris.data', header=None) # 加载Iris数据集作为DataFrame对象
X = df.iloc[:, [0, 2]].values # 取出2个特征，并把它们用Numpy数组表示
plt.scatter(X[:50, 0], X[:50, 1],color='red', marker='o', label='setosa') # 前50个样本的散点图
plt.scatter(X[50:100, 0], X[50:100, 1],color='blue', marker='x', label='versicolor') # 中间50个样本的散点图
plt.scatter(X[100:, 0], X[100:, 1],color='green', marker='+', label='Virginica') # 后50个样本的散点图
plt.xlabel('petal length') #花瓣长度
plt.ylabel('sepal length') #花萼长度
plt.legend(loc=2) # 把说明放在左上角，具体请参考官方文档
plt.show()

从上图我们可以看出，数据集是线性可分的。为了更好地演示效果，我只用setosa和versicolor这两个类别。
接下来，我们要用强大的scikit-learn来拟合模型。

# -*- coding: utf-8 -*-
"""
a test
"""
from sklearn import datasets
import numpy as np
from sklearn.cross_validation import train_test_split
iris = datasets.load_iris() # 由于Iris是很有名的数据集，scikit-learn已经原生自带了。
X = iris.data[:, [2, 3]]
y = iris.target # 标签已经转换成0，1，2了
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=0) # 为了看模型在没有见过数据集上的表现，随机拿出数据集中30%的部分做测试
# 为了追求机器学习和最优化算法的最佳性能，我们将特征缩放
from sklearn.preprocessing import StandardScaler
sc = StandardScaler()
sc.fit(X_train) # 估算每个特征的平均值和标准差
sc.mean_ # 查看特征的平均值，由于Iris我们只用了两个特征，所以结果是array([ 3.82857143,  1.22666667])
sc.scale_ # 查看特征的标准差，这个结果是array([ 1.79595918,  0.77769705])
X_train_std = sc.transform(X_train)
# 注意：这里我们要用同样的参数来标准化测试集，使得测试集和训练集之间有可比性
X_test_std = sc.transform(X_test)
# 训练感知机模型
from sklearn.linear_model import Perceptron
# n_iter：可以理解成梯度下降中迭代的次数
# eta0：可以理解成梯度下降中的学习率
# random_state：设置随机种子的，为了每次迭代都有相同的训练集顺序
ppn = Perceptron(n_iter=40, eta0=0.1, random_state=0)
ppn.fit(X_train_std, y_train)
# 分类测试集，这将返回一个测试结果的数组
y_pred = ppn.predict(X_test_std)
# 计算模型在测试集上的准确性，我的结果为0.9，还不错
from sklearn import metrics
print metrics.accuracy_score(y_test, y_pred)

scikit-learn实现逻辑回归

# -*- coding: utf-8 -*-
"""
scikit-learn实现逻辑回归
"""
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.colors import ListedColormap
def plot_decision_regions(X, y, classifier, test_idx=None, resolution=0.02):
    # setup marker generator and color map
    markers = ('s', 'x', 'o', '^', 'v')
    colors = ('red', 'blue', 'lightgreen', 'gray', 'cyan')
    cmap = ListedColormap(colors[:len(np.unique(y))])
    # plot the decision surface
    x1_min, x1_max = X[:, 0].min() - 1, X[:, 0].max() + 1
    x2_min, x2_max = X[:, 1].min() - 1, X[:, 1].max() + 1
    xx1, xx2 = np.meshgrid(np.arange(x1_min, x1_max, resolution), np.arange(x2_min, x2_max, resolution))
    Z = classifier.predict(np.array([xx1.ravel(), xx2.ravel()]).T)
    Z = Z.reshape(xx1.shape)
    plt.contourf(xx1, xx2, Z, alpha=0.4, cmap=cmap)
    plt.xlim(xx1.min(), xx1.max())
    plt.ylim(xx2.min(), xx2.max())
    # plot class samples
    for idx, cl in enumerate(np.unique(y)):
        plt.scatter(x=X[y == cl, 0], y=X[y == cl, 1],alpha=0.8, c=cmap(idx),marker=markers[idx], label=cl)
    # highlight test samples
    if test_idx:
        X_test, y_test = X[test_idx, :], y[test_idx]
        plt.scatter(X_test[:, 0], X_test[:, 1], c='', alpha=1.0, linewidth=1, marker='o', s=55, label='test set')
from sklearn import datasets
import numpy as np
from sklearn.cross_validation import train_test_split
iris = datasets.load_iris()
X = iris.data[:, [2, 3]]
y = iris.target
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.3, random_state=0)
from sklearn.preprocessing import StandardScaler
sc = StandardScaler()
sc.fit(X_train)
X_train_std = sc.transform(X_train)
X_test_std = sc.transform(X_test)
X_combined_std = np.vstack((X_train_std, X_test_std))
y_combined = np.hstack((y_train, y_test))
from sklearn.linear_model import LogisticRegression
lr = LogisticRegression(C=1000.0, random_state=0)
lr.fit(X_train_std, y_train)
lr.predict_proba(X_test_std[0,:]) # 查看第一个测试样本属于各个类别的概率
plot_decision_regions(X_combined_std, y_combined, classifier=lr, test_idx=range(105,150))
plt.xlabel('petal length [standardized]')
plt.ylabel('petal width [standardized]')
plt.legend(loc='upper left')
plt.show()