Question

什么是回归分析？回归分析是一项预测性的建模技术。它的目的是通过建立模型来研究因变量和自变量之间的显著关系，即多个自变量对（一个）因变量的影响强度，预测数值型的目标值。

回归分析在管理、经济、社会学、医学、生物学等领域得到了广泛的应用，这种技术的起源最早可以追溯到达尔文（Charles Darwin）时期。当时，他的表兄弟Francis Galton正致力于研究父代豌豆种子尺寸对子代豌豆尺寸的影响，采用了回归分析，并同时在多项研究中都注意到这个现象。在19世纪高斯系统地提出最小二乘估计后，回归分析开始蓬勃发展。目前，回归分析的研究范围可分为几大板块，如下：

常用的回归分析技术是线性回归、逻辑回归、多项式回归和岭回归等。作为入门书籍，在此主要介绍前两种模型的原理和具体实现。

7.1.1 线性回归

线性回归是最简单的回归模型，如图7-1所示。它的目的是：在自变量（输入数据）仅一维的情况下，找出一条最能够代表所有观测样本的直线（估计的回归方程）。在自变量（输入数据）高于一维的情况下，找到一个超平面使得数据点距离这个平面的误差（Residuals）最小。而前者的情况在高中数学课本就已经学过，它的解法是普通最小二乘法（Ordinary Least Squares，OLS）。

图7-1 线性回归示意图

线性回归模型的公式如下所示。

f（x）=w₁ x₁ +w₂ x₂ +…+w_n x_n +b×1

其中，权重W_i 和常数项b是待确定的。这意味着将输入的自变量按一定比例加权求和，得到预测值输出。

1.算法推导

重述一遍普通最小二乘法是不必要的，我们将详细介绍更容易推广到高维情况的矩阵推导过程。一般来说，数据挖掘算法都会涉及矩阵的表示、运算和结果推导。这是因为矩阵的本质是一张数表，与常见的数据格式很贴合；同时，利用抽象的矩阵来表示算法推导过程非常简洁。读者将会从后续几个算法的推导中逐步感受到。

假设矩阵X_m×n 代表m个样本n维特征的数据对应的矩阵，且X_m×n 的列向量线性无关。通常，我们会在X_m×n 的最后一列添加一列全为1的向量，以对应上述公式中提到的截距。此时，权重向量ω_（n+1）×1 =（w₁ ，w₂ ，…，w_n ，b）^T ，目标是计算权重向量使得预测值Xω与真实值y的均方误差最小。公式如下：

min E（ω）=|Xω-y|²

其中，E（ω）可化简为如下公式：

E（ω）=（Xω-y）^T （Xω-y）=ω^T X^T Xω-2y^T Xω+y^T y

由于E（ω）是一个凸函数，因此其理论最小值出现在偏导数全为0处。注意，此处是对向量求导，要求一些矩阵分析或矩阵微积分的计算经验。

因为X_m×（n+1） 的列向量线性无关，所以X^T X总是可逆的。化简上式求得：

ω=（X^T X）^-1 X^T y

显然，如果X_m×n 退化为一个列向量 ，整个推导过程与简单线性回归f（x）=wx+b×1对应的普通最小二乘法无异。

2.算法实现

通过前面的学习我们知道，Python有强大的第三方扩展模块sklearn（读作scikit-learn），实现了绝大部分的数据挖掘基础算法，包括线性回归。下面我们将通过举例说明，如何使用sklearn快速实现线性回归模型。这个例子是经典的波士顿房价预测问题^[1] ，如表7-1所示。

表7-1 波士顿房价预测的部分数据

*数据详见：示例程序/data/BostonHousePricing.csv

直觉告诉我们：数据表7-1中住宅平均房间数列与最终房价一般是成正比的，具有某种线性关系。我们利用线性回归来验证想法。同时，作为一个二维的例子，可以在平面上绘出图形，进一步观察图形，如代码清单7-1所示，效果见图7-2。

代码清单7-1 线性回归sklearn实现

# -*- coding:utf-8 -*-
import numpy as np
import matplotlib.pyplot as plt
from sklearn.datasets import load_boston
from sklearn.linear_model import LinearRegression
boston = load_boston()
print boston.keys()
# result:
# ['data', 'feature_names', 'DESCR', 'target']
print boston.feature_names
# result:
# ['CRIM' 'ZN' 'INDUS' 'CHAS' 'NOX' 'RM' 'AGE' 'DIS' 'RAD' 'TAX' 'PTRATIO' 'B' 'LSTAT']
# print boston.DESCR          # 取消注释并运行，可查看数据说明文档
x = boston.data[:, np.newaxis, 5]
y = boston.target
lm = LinearRegression()          # 声明并初始化一个线性回归模型的对象
lm.fit(x, y)               # 拟合模型，或称为训练模型
print u'方程的确定性系数
(R^2): %.2f' % lm.score(x, y)
# result: 方程的确定性系数
(R^2): 0.48
plt.scatter(x, y, color='green')     # 显示数据点
plt.plot(x, lm.predict(x), color='blue', linewidth=3)     # 画出回归直线
plt.xlabel('Average Number of Rooms per Dwelling (RM)')
plt.ylabel('Housing Price')
plt.title('2D Demo of Linear Regression')
plt.show()

*代码详见：示例程序/code/7-1-1.py

7.1.2 逻辑回归

分类和回归二者不存在不可逾越的鸿沟。以波士顿房价预测为例：如果将房价按高低分为“高级”、“中级”和“普通”三个档次，那么这个预测问题也属于分类问题。当然，我们要注意保持训练集和测试集的一致性。如果是住房档次的类别预测问题，我们首先应该将原始数据按不同价格区间分档，改写数据后再进行后续步骤。

图7-2 线性回归示例图

准确地说，逻辑回归 （Logistic Regression）是对数几率回归，属于广义线性模型 （GLM），它的因变量一般只有0或1。读者需要明确一件事情：线性回归并没有对数据的分布进行任何假设，而逻辑回归隐含了一个基本假设h：每个样本均独立服从于伯努利分布 （0-1分布）。伯努利分布属于指数分布族，这个大家族还包括：高斯（正态）分布、多项式分布、泊松分布、伽马分布、Dirichlet分布等。

事实上，假设数据服从某类指数分布，我们可以由线性模型拓展出一类广义线性模型，即通过非线性的关联函数 （Link Function）将线性函数映射到其他空间上，从而大大扩大了线性模型本身可解决的问题范围。根据数理统计的基础知识，指数分布的概率密度函数可抽象地表示为：

其中，η是待定的参数，T（y）是充分统计量，通常T（y）=y。

而伯努利分布的概率密度函数为：

其中，h是由假设衍生的关联函数，y的可能取值为0或1。值得注意的是，当y=1时，p（y；h）=h。也就是说，h表征样本属于正类（类别“1”）的概率。对照如下公式，令 ，可得：

这便是大名鼎鼎的Logistic函数，亦称Sigmoid函数，如图7-3所示。因为它的函数形如字母“S”。

图7-3 Logistic函数

观察图像可知，当指数分布的自然参数η在（-∞，+∞）变化时，h的取值范围恰好为（0，1）。由于h表征样本属于正类（类别“1”）的概率，通常将h大于某个阈值（如0.5）的样本预测为“属于正类（1）”，否则预测结果为“属于负类（0）”。

由基本假设，η=w^T x。给定x，目标函数为

上式表明，我们最终目标是根据数据确定合适的一组权重w。在对原始输入进行加权组合之后，通过关联函数做非线性变换，得到的结果表示样本x属于正类的概率。因此，关联函数亦称为“激活函数”，如同神经元接受到足够的刺激，才会变得兴奋。

通常，确定权重w采用的方法是：最大似然估计（Maximum Likelihood Estimation）。这在任何一本数理统计教材中都有讲到。在此，通过一个简单的例子来重温它的思想。

例如一枚来自赌场的硬币，有“正”“反”两面。这枚硬币对应了一个未知的参数：η₀ ，即抛出正面的概率。于是你重复抛这枚硬币10000次，其中有8000次都是正面。基于这个事实（数据），我们不会再愿意相信“硬币抛出正反两面的概率都是0.5”这条假设，而更倾向于认为“这枚硬币抛出正面的概率是0.8”，即η₀ =0.8。这便是最大似然估计的思想：基于事实、样本结果或数据等，来估计（确认）一个概率模型的关键参数。这未必是完全正确的，却是迄今为止最好的。

算法实现

工程中求解逻辑回归更倾向于选择一些迭代改进的算法，如牛顿方法、梯度下降等。它们会直接对解空间进行部分搜索，找到合适的结果便停止寻优。建议读者在入门时首先掌握scikit-learn中的逻辑回归实现算法，如代码清单7-2所示。

代码清单7-2 逻辑回归sklearn实现

# -*- coding:utf8 -*-
import pandas as pd
from sklearn.linear_model import LogisticRegression, RandomizedLogisticRegression
from sklearn.cross_validation import train_test_split
# 导入数据并观察
data = pd.read_csv('../data/LogisticRegression.csv', encoding='utf-8')
# print data.head(5)     # 查看数据框的头五行
# 将类别型变量进行独热编码
 (one-hot encoding)
data_dum = pd.get_dummies(data, prefix='rank', columns=['rank'], drop_first=True)
print data_dum.tail(5)     # 查看数据框的最后五行
# result:
#     admit  gre   gpa  rank_2  rank_3  rank_4
# 395      0  620  4.00     1.0     0.0     0.0
# 396      0  560  3.04     0.0     1.0     0.0
# 397      0  460  2.63     1.0     0.0     0.0
# 398      0  700  3.65     1.0     0.0     0.0
# 399      0  600  3.89     0.0     1.0     0.0
# 切分训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(data_dum.ix[:, 1:], data_dum.ix[:, 0], test_size=.1, random_state=520)
lr = LogisticRegression()     # 建立
LR模型
lr.fit(X_train, y_train)     # 用处理好的数据训练模型
print '逻辑回归的准确率为：
{0:.2f}%'.format(lr.score(X_test, y_test) *100)

*代码详见：示例程序/code/7-1-2.py

[1] 在sklearn的安装文件中，已包含此问题对应的数据集。