Logistic Regression原理及Python实现

1. 问题引入

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

2. Logistic Regression

继续我们上面的预测患者是否会犯心脏病的问题。为了最终的预测，我们可以对患者的身体状况x=(x0,x1,x2⋯,xd)按照不同的权重打分：

s=∑i=0dwixi

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

θ(s)=es1+es=11+e−s

也就是说，对数几率回归通过h(x)=11+exp(−wTx)来逼近目标函数f(x)=P(+1|x)，显然有：

P(y=1|x)=f(x)P(y=0|x)=1−f(x)

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

假设现在有一个数据集D={(x1,y1),(x2,y2),⋯,(xN,yn)}，生成该数据集的概率为：

∏i=1N[P(yi|xi)]yi[1−P(yi|xi)]1−yi

由于P(yi|xi)未知，所以我们用h(xi)来代替它：

∏i=1N[h(xi)]yi[1−h(xi)](1−yi)

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

L(w)=∑n=1N[yilogh(xi)+(1−yi)log(1−h(xi))]=∑n=1N[yilogh(xi)1−h(xi)+log(1−h(xi))]=∑n=1N[yilogexp(w⋅xi)/(1+exp(w⋅xi))1/(1+exp(w⋅xi))+log(1−h(xi))]=∑n=1N[yi(w⋅xi)+log(1+exp(w⋅xi))]

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

∇L(w)=∑n=1N[yi⋅xi−11+exp(w⋅xi)⋅exp(w⋅xi)⋅xi]=∑n=1N[yi⋅xi−exp(w⋅xi)1+exp(w⋅xi)⋅xi]=∑n=1N[xi⋅(yi−exp(w⋅xi)1+exp(w⋅xi))]=∑n=1N[xi⋅(yi−h(w⋅xi))]

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

w^(k+1)←w^(k)+η⋅∇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×3的矩阵，其中第一列全是1，表示模型中的常数（类似于线性回归中的b），而label_mat是一个100×1的矩阵。下面我们先实现sigmoid函数：

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

接下来我们实现对数几率回归。logistic_regression函数一共有四个参数：第一个参数代表输入的数据；第二个参数表示输入数据的标签；参数alpha代表学习率，默认为0.001；最后一个参数代表梯度上升的迭代次数，默认为500次。我们首先初始化变量weight为一个3×1的全1矩阵，然后在每一轮迭代中计算h(x⋅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()

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