回归模型原理与应用 （波士顿房价预测）

回归一词由弗朗西斯·高尔顿爵士（1822-1911）提出，他发现父母一高一矮的人，身高区域父母身高之间，这种现象被他称为“向均值回归”。回归是研究自变量X和因变量Y之间的关系。X与Y之间的关系可以用回归函数表示，所以回归问题的估计可以视为函数的拟合。本问假设X与Y是线性关系，将为读者介绍线性回归和logistic回归，详细讲解最小二乘法，以及结合实际问题进行应用。

1.1 理论模型

比如说降雨量和粮食产量的关系，房子面积和租金的关系等等，这是X只有一个特征（如降雨量和房子面积）的情况。再比如地段、房子面积、房子类型、银行贷款利率等于房价的关系，这是一个多维特征X与Y的关系。具体模型如下：
给定随机样本$(x_i,y_i)，i=1,...,m$，$x_i$有n个属性。线性模型为：$$\begin{gathered} Y_i={\beta }_0+{\beta }_1{x}_{i1}+{\beta }_2{x}_{i2}+\cdots +{\beta }_n{x}_{in}+{ϵ}_i,i=1,2,\dots ,m \\ \mathbb{E}({ϵ}_i|x_i)=0，\mathbb{V}({ϵ}_i|x_i)={\sigma }^2 \end{gathered}$$
这里的${\beta }_0，{\beta }_1，\dots ,{\beta }_n$是模型中需要估计的参数，$ϵ$是残差项。这里我们需要要求残差项服从正态分布。

思考：为什么会有这个残差项呢？
因为我们假设除了X的影响，还会有其他的因素影响Y的结果。就以之前降雨量和粮食产量的例子，根据生活常识，除了降雨量之外，肯定还有其他未知的因素影响到粮食产量，并且这个因素对粮食产量的影响是随机的，这个未知因素的影响就了残差项中。此外，对于不随机的因素，我们用参数${\beta }_0$表示偏移。就比如二维情况下的截距。

因为有n维的特征，按照上述的表示显得比较繁琐，我们不妨将其向量化，显得更加清晰明了。
记:
$$Y={\left(y_1,y_2,y_3,\dots ,y_m\right)}^T$$$$\beta ={\left({\beta }_0,{\beta }_1,{\beta }_2,\dots ,{\beta }_m\right)}^T$$$$ϵ={\left({ϵ}_1,{ϵ}_2,{ϵ}_3,\dots ,{ϵ}_m\right)}^T$$$$X=\left(\begin{array}{cccccc}1& x_{11}& x_{12}& x_{13}& \dots & x_{1n}\\ 1& x_{21}& x_{22}& x_{23}& \dots & x_{2n}\\ ⋮& ⋮& ⋮& ⋮& \dots & ⋮\\ 1& x_{m1}& x_{m2}& x_{m3}& \dots & x_{mn}\end{array}\right)$$
则线性模型可以表示为：
$$Y=X\beta +ϵ$$
这里需要注意，因为参数β中含有${\beta }_0$，所以需要在X中添加全为1的列。

1.2 数据和估计

在使用样本数据集估计上述模型中的参数时，我们首先要问的一个问题是，对于数据有什么要求？对机器学习有一些了解得读者肯定知道数据量太小会造成过拟合。那么怎么才算数据量太少呢？假设我们不考虑残差项，那么线性模型就可以变成成一个未知数为β的线性方程组求解问题，m个样本代表着m个方程组。在线性代数中，我们学习过方程组解与方程组秩的关系。考虑到我们是随机采样的，所以样本数可以近似看做秩。那么，当样本数小于参数个数时，方程就有无数多个解，这样就造成了过拟合。当样本数大于参数个数时，因为我们很多时候是很难把所有因素考虑进去，或者包含随机因素，系数矩阵是相容的。通过引入残差项，最下化残差来近似估计参数。

思考
形如下式的模型是否为线性模型?
$$Y_i={\beta }_0+{\beta }_1{x}_{i1}^1+{\beta }_2{x}_{i2}^2+\cdots +{\beta }_n{x}_{in}^n+{ϵ}_i,i=1,2,\dots ,m$$是线性模型，这里所说的线性不需要是自变量的线性函数。线性这里指的$Y_i$关于参数β是线性的。如上面近似看做线性方程组的理解方法，$x_i$的非线性只改变系数矩阵，要求解的仍然是个线性方程组。

因为引入了残差，所以不可以用求解线性方程组的方法准确计算参数。这里我们就要采用最小二乘法对参数进行估计。
所谓最小二乘，就是最小化残差的平方，换句话说就是最小化理论值和估计值的平方和。即：
$$\begin{array}{rl}\beta & =\underset{\beta }{\mathrm{argmin}}\sum _{i=1}^m(Y_i-X_i\beta {)}^2\\ & =\underset{\beta }{\mathrm{argmin}}\sum _{i=1}^m(Y-X\beta {)}^T(Y-X\beta )\end{array}$$
这是一个二次优化问题，可以比较容易求解得出：
当$X^{T}X$是可逆的时：
$$\hat{\beta }=(X^{T}X{)}^{-1}{X}^{T}Y$$

1.3 模型选择

对于线性模型，存在一个模型参数选择的问题。比如说：关于房价预测问题，我们可能会找到很多可能相关的特征，如：距市中心距离、房间大小、修筑时间、当前股市、城市人口等等诸多因素。但是很多因素可能关联十分的低，那么我们要不要去除掉这些因素呢？去掉会更好一些，一方面特征数少不容易发生过拟合，预测效果更好，另一方面模型简单。那么我们应该遵从什么方式选择特征呢？可以想到打分的方式，选择得分高的模型。下面介绍三种打分都是遵循：
交叉验证
交叉验证的方法是最小化下式：
$$R(s)=\sum _i^m(Y_i-{\hat{Y}}_i)$$s是样本集的子集，其$R(s)$表示以s为训练集的风险估计。等式的右边表示把第i个样本删去后估计的参数所预测的结果和实际结果误差的平方和。这种方法在机器学习中叫做留一法。这一般需要很大的计算量，我们有时会采用k折法。即将样本集分为k等分，然后求把第i组样本删去后估计的参数所预测的结果和实际结果误差的平方和。

AIC(Akaika Information Critirion)法
AIC法是要使得下式达到最大：
$$AIC(s)=l_s-|s|$$这里我们设r是选取一部分特征的X的子集，$|s|$表示要估计参数个数，也就是训练集的维度。$l_s$表示根据r极大似然估计值。对于线性回归的极大似然值近似等于最小二乘估计的方差和。在下一小节会对此进行介绍。

BIC(Bayes Information Critirion)法
BIC法是要使得下式达到最大：
$$BIC(s)=l_s-\frac{|s|}{2}\log m$$上上式中的m指的是训练样本的数目。BIC模型与AIC模型的思路一致，不过加强了对复杂度的惩罚。BIC法有一种贝叶斯解释，还有一种信息论的最小长度的解释。下面介绍贝叶斯解释，假设模型的集合 S = { s 1 , s 2 , … , s p } S=\{ s_{1},s_{2},\dots,s_{p} \} S={s1​,s2​,…,sp​}，假设每个模型的先验概率都为$\frac{1}{p}$。可以证明模型的后验概率为：
$$\mathbb{P}(s_j|数据集)=\frac{e^{BIC(s_j)}}{{\sum }_r^p{e}^{BIC(s_r)}}$$因此我们选择使得后验概率最大的模型。更多的关于BIC的知识详见[1]。
上述的三种模型选择的方法，其中交叉验证的使用率最高，因为当模型十分复杂的时候，后两种方法不太实用。后两种方法中，BIC相对AIC，对复杂度的惩罚力度更大。在实际应用中，参数很多的情况下，我们要计算每一种特征排列组和的结果是十分耗费计算资源的。一般我们不遍历每一种特征的组合，而是采用前向或者(后向)逐步回归的方法。即从只有使得打分最高1个特征，选择加一个特征是打分变得最高的加入特征组合，一直递归下去知道模型得分不再提高。

1.3 基于波士顿房价预测的实战演练

实战演练可以让我们对于知识的理解更加深入。本节我们使用波士顿房价预测的数据集对线性回归模型进行测试。

1. 导入需要的包

import numpy as np           # 数学运算模块
import matplotlib.pyplot as plt # 可视化
from sklearn.datasets import load_boston # 数据集
from sklearn.linear_model import LinearRegression # 线性回归模块

2. 导入数据集

boston = load_boston()
print(boston.DESCR)

运行结果

---

Data Set Characteristics:

:Number of Instances: 506 

:Number of Attributes: 13 numeric/categorical predictive. Median Value (attribute 14) is usually the target.

:Attribute Information (in order):
    - CRIM     per capita crime rate by town
    - ZN       proportion of residential land zoned for lots over 25,000 sq.ft.
    - INDUS    proportion of non-retail business acres per town
    - CHAS     Charles River dummy variable (= 1 if tract bounds river; 0 otherwise)
    - NOX      nitric oxides concentration (parts per 10 million)
    - RM       average number of rooms per dwelling
    - AGE      proportion of owner-occupied units built prior to 1940
    - DIS      weighted distances to five Boston employment centres
    - RAD      index of accessibility to radial highways
    - TAX      full-value property-tax rate per $10,000
    - PTRATIO  pupil-teacher ratio by town
    - B        1000(Bk - 0.63)^2 where Bk is the proportion of blacks by town
    - LSTAT    % lower status of the population
    - MEDV     Median value of owner-occupied homes in $1000's

:Missing Attribute Values: None

:Creator: Harrison, D. and Rubinfeld, D.L.

This is a copy of UCI ML housing dataset.
https://archive.ics.uci.edu/ml/machine-learning-databases/housing/

This dataset was taken from the StatLib library which is maintained at Carnegie Mellon University.

The Boston house-price data of Harrison, D. and Rubinfeld, D.L. ‘Hedonic
prices and the demand for clean air’, J. Environ. Economics & Management,
vol.5, 81-102, 1978. Used in Belsley, Kuh & Welsch, ‘Regression diagnostics
…’, Wiley, 1980. N.B. Various transformations are used in the table on
pages 244-261 of the latter.

The Boston house-price data has been used in many machine learning papers that address regression
problems.

… topic:: References

• Belsley, Kuh & Welsch, ‘Regression diagnostics: Identifying Influential Data and Sources of Collinearity’, Wiley, 1980. 244-261.
• Quinlan,R. (1993). Combining Instance-Based and Model-Based Learning. In Proceedings on the Tenth International Conference of Machine Learning, 236-243, University of Massachusetts, Amherst. Morgan Kaufmann.

3. 初步认识数据集

feature_names = list(boston.feature_names)
X = boston.data
Y = boston.target
print(boston.feature_names)
print("特征个数:",X.shape[1])
print("样本数:",Y.shape[0])
# print(boston)
print("价格最小值: $1000's{:,.2f}".format(np.min(Y)))
print("价格最大值: $1000's{:,.2f}".format(np.max(Y)))
print("平均价格: $1000's{:,.2f}".format(np.mean(Y)))
print("价格中位数: $1000's{:,.2f}".format(np.median(Y)))
print("价格标准差: $1000's{:,.2f}".format(np.std(Y)))

运行结果：

[‘CRIM’ ‘ZN’ ‘INDUS’ ‘CHAS’ ‘NOX’ ‘RM’ ‘AGE’ ‘DIS’ ‘RAD’ ‘TAX’ ‘PTRATIO’ ‘B’ ‘LSTAT’]
特征个数: 13
样本数: 506
价格最小值: $1000’s5.00
价格最大值: $1000’s50.00
平均价格: $1000’s22.53
价格中位数: $1000’s21.20
价格标准差: $1000’s9.19

4. 计算模型打分
   这里采用了10折交叉验证的方法对模型进行打分。

def get_score(X_s):
    cv = KFold(n_splits=10, random_state=0)
    # print(sorted(sm.SCORERS.keys()))
    lm = LinearRegression()
    scoring = ['neg_mean_squared_error', 'neg_mean_absolute_error'] # precision_macro为精度，recall_macro为召回率
    scores = cross_val_score(lm, X_s, Y, scoring="neg_mean_squared_error",cv=cv)
    return scores.mean()
print("10折交叉验证平均均方误差：",-get_score(X))

10折交叉验证平均均方误差：34.705255944524964

5.前向递归选择特征
方法描述见模型选择

all_features = feature_names.copy()
all_score = -1000
features_index = []
def choose_features(all_score):
    global features_index
    temp_list = []
    scores = []
    for i,feature in enumerate(all_features):
        temp_features = features_index + [feature_names.index(feature)]
        score = get_score(X[:,temp_features])
        scores.append(score)
        temp_list.append(temp_features)
    if (max(scores) - all_score) > 0.001:
        features_index = temp_list[scores.index(max(scores))]
        all_features.pop(scores.index(max(scores)))
        all_score = max(scores)
        print("10折交叉验证平均均方误差：",-max(scores))
        return choose_features(all_score)
    else:
        return features_index
choose_features(all_score)   
print("特征选择：",boston.feature_names[features_index])

10折交叉验证平均均方误差： 41.828958072164035
10折交叉验证平均均方误差： 35.24329476929434
10折交叉验证平均均方误差： 33.160249867250386
10折交叉验证平均均方误差： 32.51383989910426
10折交叉验证平均均方误差： 31.91687506868982
10折交叉验证平均均方误差： 31.609745485338117
10折交叉验证平均均方误差： 31.39917693777587
10折交叉验证平均均方误差： 31.160306011122884
10折交叉验证平均均方误差： 31.12986419207171
特征选择： [‘LSTAT’ ‘PTRATIO’ ‘RM’ ‘DIS’ ‘NOX’ ‘CHAS’ ‘CRIM’ ‘ZN’ ‘INDUS’]

6. 结果分析
   通过基于交叉验证的前向递归算法，我们选取了10个特征，但从递归过程来看，主要起到作用的是前四个特征，分别是：
   ‘LSTAT’：该地区有多少百分比的业主属于是低收入阶层；
   ‘PTRATIO’ ：该地区的中学和小学里，学生和老师的数目比；
   ‘RM’：该地区中每个房屋的平均房间数量；
   ‘DIS’：到五个波士顿就业中心的加权距离。
   看一下单独与每个特征与房价的散点图。

# 使输出的图像以更高清的方式显示
%config InlineBackend.figure_format = 'retina'
# 调整图像的宽高
plt.figure(figsize=(16, 4))
for i, key in enumerate([12, 10, 5, 7,]):
    plt.subplot(1, 4, i+1)
    plt.xlabel(feature_names[key])
    plt.scatter(X[:,key], Y, alpha=0.5)

运行结果