Logistic Regression原理及Python实现

1. 问题引入

相信大家都接触过分类问题，尤其是二元分类。例如现在有一些患者（训练集）的身体情况以及是否患有心脏病的数据，要求我们根据这些数据来预测其他患者（测试集）是否患有心脏病。这是比较简单的一个二元分类问题，使用线性分类器或许会取得不错的效果。但在实际生活中，我们感兴趣的往往不是其他患者是否会犯病，而是他犯心脏病的概率是多少。很直观的想法是收集患者犯病的概率，然后利用回归模型进行概率预测。但是我们并不能直接收集到患者犯病的概率，只能知到患者最后到底犯没犯病。也就是说，我们输入的标签是离散的类别型数据，但是期望得到数值型的概率。在机器学习中，这叫做“软性二元分类”（Soft Binary Classification），这类问题与普通的分类问题所要的数据相同，但是会得到不同的目标函数。接下来，我们就介绍一种可以解决“软性二元分类”问题的算法，那就是Logistic Regression（后文译为“对数几率回归）。

2. Logistic Regression

继续我们上面的预测患者是否会犯心脏病的问题。为了最终的预测，我们可以对患者的身体状况$\mathbf{x}=(x_0,x_1,x_2\cdots ,x_d)$按照不同的权重打分：

$$s=\sum _{i=0}^d{w}_i{x}_i$$

然后用对数几率函数$\theta (s)$将分数转化成一个概率估计。对数几率函数的形式为：

$$\theta (s)=\frac{e^s}{1+e^s}=\frac{1}{1+e^{-s}}$$

也就是说，对数几率回归通过$h(\mathbf{x})=\frac{1}{1+\exp (-{\mathbf{w}}^T\mathbf{x})}$来逼近目标函数$f(\mathbf{x})=P(+1|\mathbf{x})$，显然有：

$$\begin{gathered} P(y=1|\mathbf{x})=f(\mathbf{x}) \\ P(y=0|\mathbf{x})=1-f(\mathbf{x}) \end{gathered}$$

现在我们要面临的问题是如何对$\mathbf{w}$进行估计。在线性回归模型中，我们可以通过最小化均方误差来估计回归系数。但是在这里，好像并没有什么可以衡量误差的东西，所以我们要另辟蹊径。下面我们使用极大似然估计（Maximum Likelihood Estimation，简称MLE）来估计最优的回归系数，它是根据数据采样来估计概率分布参数的经典方法。

假设现在有一个数据集D={(x1,y1),(x2,y2),⋯ ,(xN,yn)}D = \{ (\mathbf{x}_1, y_1),(\mathbf{x}_2,y_2),\cdots,(\mathbf{x}_N,y_n) \}D={(x1​,y1​),(x2​,y2​),⋯,(xN​,yn​)}，生成该数据集的概率为：

$$\prod _{i=1}^N[P(y_i|{\mathbf{x}}_i){]}^{y_i}[1-P(y_i|{\mathbf{x}}_i){]}^{1-y_i}$$

由于$P(y_i|{\mathbf{x}}_i)$未知，所以我们用$h({\mathbf{x}}_i)$来代替它：

$$\prod _{i=1}^N[h({\mathbf{x}}_i){]}^{y_i}[1-h({\mathbf{x}}_i){]}^{(1-y_i)}$$

由于连乘不好优化，所以要将上式变成连加的形式，对上式取对数得对数似然函数：

$$\begin{array}{rl}L(w)& =\sum _{n=1}^N\left[y_i\log h(x_i)+(1-y_i)\log (1-h(x_i))\right]\\ & =\sum _{n=1}^N\left[y_i\log \frac{h(x_i)}{1-h(x_i)}+\log (1-h(x_i))\right]\\ & =\sum _{n=1}^N\left[y_i\log \frac{\exp (w\cdot x_i)/(1+\exp (w\cdot x_i))}{1/(1+\exp (w\cdot x_i))}+\log (1-h(x_i))\right]\\ & =\sum _{n=1}^N\left[y_i(w\cdot x_i)+\log (1+\exp (w\cdot x_i))\right]\end{array}$$

对$L(w)$求极大值即可得到$w$的估计值。这样，问题就变成了以对数似然函数为目标函数的最优化问题。求极大值可以采用梯度上升法。首先求得$L(w)$在$w$上的梯度：

$$\begin{array}{rl}\nabla L(w)& =\sum _{n=1}^N\left[y_i\cdot x_i-\frac{1}{1+\exp (w\cdot x_i)}\cdot \exp (w\cdot x_i)\cdot x_i\right]\\ & =\sum _{n=1}^N\left[y_i\cdot x_i-\frac{\exp (w\cdot x_i)}{1+\exp (w\cdot x_i)}\cdot x_i\right]\\ & =\sum _{n=1}^N\left[x_i\cdot \left(y_i-\frac{\exp (w\cdot x_i)}{1+\exp (w\cdot x_i)}\right)\right]\\ & =\sum _{n=1}^N\left[x_i\cdot \left(y_i-h(w\cdot x_i)\right)\right]\end{array}$$

那么，第k轮迭代时$w$的极大似然估计值为：

$${\hat{w}}^{(k+1)}\leftarrow {\hat{w}}^{(k)}+\eta \cdot \nabla L(w^{(k)})$$

经过多轮迭代后，我们便可以求得$w$的极大似然估计值。下面我们使用Python来实现对数几率回归这一算法。

首先，我们引入需要用到的模块，并加载数据集。

import numpy as np
import pandas as pd
from pandas import Series, DataFrame
import matplotlib.pyplot as plt
%matplotlib inline

df = pd.read_table('testSet.txt',header=None)
df.shape

(100, 3)

df.head()

	0	1	2
0	-0.017612	14.053064	0
1	-1.395634	4.662541	1
2	-0.752157	6.538620	0
3	-1.322371	7.152853	0
4	0.423363	11.054677	0

data_mat = DataFrame({0:1, 1:df[0], 2:df[1]})
data_mat.head()

	0	1	2
0	1	-0.017612	14.053064
1	1	-1.395634	4.662541
2	1	-0.752157	6.538620
3	1	-1.322371	7.152853
4	1	0.423363	11.054677

label_mat = DataFrame(df[2])
label_mat.head()

	2
0	0
1	1
2	0
3	0
4	0

通过观察数据集，我们可以看出data_mat是一个$100\times 3$的矩阵，其中第一列全是$1$，表示模型中的常数（类似于线性回归中的$b$），而label_mat是一个$100\times 1$的矩阵。下面我们先实现sigmoid函数：

def sigmoid(x_arr):
    return 1.0/(1+np.exp(-x_arr))

接下来我们实现对数几率回归。logistic_regression函数一共有四个参数：第一个参数代表输入的数据；第二个参数表示输入数据的标签；参数alpha代表学习率，默认为0.001；最后一个参数代表梯度上升的迭代次数，默认为500次。我们首先初始化变量weight为一个$3\times 1$的全1矩阵，然后在每一轮迭代中计算$h(x\cdot w)$的值，然后计算h与y_mat的差值error，error与x_mat相乘就是梯度，梯度与alpha相乘之后就可以更新weight。注意这里全都是矩阵运算。

def logistic_regression(x_arr, y_arr, alpha=.001, max_iter=500):
    x_mat = np.mat(x_arr)
    y_mat = np.mat(y_arr)
    weight = np.ones((x_mat.shape[1],1))
    for k in range(max_iter):
        h = sigmoid(x_mat * weight)
        error = (y_mat - h)
        weight = weight + alpha * x_mat.T * error
    return weight

将data_mat和label_mat作为参数传递给logistic_regression函数，得到回归系数的最佳估计值：

weight = logistic_regression(data_mat, label_mat)
weight

matrix([[ 4.12414349],
        [ 0.48007329],
        [-0.6168482 ]])

最后，通过图形来观察我们的模型：

data_arr = np.array(data_mat)
n = data_arr.shape[0]

xcord1 = list()
ycord1 = list()
xcord2 = list()
ycord2 = list()

for i in range(n):
    if int(label_mat.iloc[i]) == 1:
        xcord1.append(float(data_arr[i][1]))
        ycord1.append(float(data_arr[i][2]))
    else:
        xcord2.append(data_arr[i,1])
        ycord2.append(data_arr[i,2])
fig = plt.figure(figsize=(16,12))
plt.scatter(xcord1,ycord1,s=30,c='red',marker='s')
plt.scatter(xcord2,ycord2,s=30,c='green')
x = np.arange(-3.0, 3.0, .1)
y = -(weight[0] + weight[1] * x) / weight[2]
plt.plot(x,y)
plt.show()

可以看到，这条直线基本上将红点和绿点全部区分开来，这说明我们的模型效果还不错。