python3机器学习sklearn之线性回归(一)

0. 引言

• 线性回归(Linear Regression)是利用数理统计中回归分析，来确定两种或两种以上变量间相互依赖的定量关系的一种统计分析方法
• 线性回归利用称为线性回归方程的最小平方函数对一个或多个自变量和因变量之间关系进行建模。 这种函数是一个或多个称为回归系数的模型参数的线性组合。只有一个自变量的情况称为简单回归,大于一个自变量情况的叫做多元回归。

本章介绍用线性模型处理回归问题。从简单问题开始：
* 先处理一个输入变量和一个输出变量的一元问题。

$$y = w_0 + w_1x$$

• 然后，我们介绍多元线性回归问题（multiple linear regression），线性约束由多个解释变量构成。
  

  $$y = w_0 + w_1x_1 + w_2x_2 + ..... + w_nx_n$$

• 紧接着，我们介绍多项式回归分析（polynomial regression问题），一种具有非线性关系的多元线性回归问题。
  

  $$y = w_0 + w_1x_1 + w_2x_2^2 + ..... + w_nx_n^n$$

• 最后，我们介绍如果训练模型获取目标函数最小化的参数值。在研究一个大数据集问题之前，我们先从一个小问题开始学习建立模型和学习算法。

2. 线性回归的实际用途

线性回归有很多实际的用途，分为以下两类：
1. 如果目标是预测或者映射，线性回归可以用来对观测数据集的y和X的值拟合出一个预测模型。当完成这样一个模型以后，对于一个新增的X值，在没有给定与它相配对的y的情况下，可以用这个拟合过的模型预测出一个y值。

1. 给定一个变量y和一些变量$X_1$, ⋯ , $X_p$,这些变量有可能与y相关，线性回归分析可以用来量化y与$X_j$间相关性的强度，评估出与y不相关的$X_j$ ，并识别出哪些$X_j$的子集包含了关于y的冗余信息。

3. 一元线性回归

上一章我们介绍过在监督学习问题中用训练数据来估计模型参数。训练数据由解释变量的历史观测值
和对应的响应变量构成。模型可以预测不在训练数据中的解释变量对应的响应变量的值。回归问题的
目标是预测出响应变量的连续值。本章我们将学习一些线性回归模型，后面会介绍训练数据，建模和
学习算法，以及对每个方法的效果评估。

实例一：计算匹萨的价格

首先，我们从简单的一元线性回归问题开始。
假设你想计算匹萨的价格。虽然看看菜单就知道了，不过也可以用机器学习方法建一个线性回归模
型，通过分析匹萨的直径与价格的数据的线性关系，来预测任意直径匹萨的价格。我们先用scikit
learn写出回归模型，然后我们介绍模型的用法，以及将模型应用到具体问题中。假设我们查到了部
分匹萨的直径与价格的数据，这就构成了训练数据，如下表所示：

3.1.可行性分析

• 简单而直观的方式是通过数据的可视化直接观察房屋成交价格与房屋尺寸间是否存在线性关系。

• 对于本实验的数据来说，散点图就可以很好的将其在二维平面中进行可视化表示。

我们可以用matplotlib画出图形：

%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib.font_manager import FontProperties
font = FontProperties(fname=r"c:\windows\fonts\msyh.ttc", size=10)

def runplt(size=None):
    plt.figure(figsize=size)
    plt.title('匹萨价格与直径数据',fontproperties=font)
    plt.xlabel('直径（英寸）',fontproperties=font)
    plt.ylabel('价格（美元）',fontproperties=font)
    plt.axis([0, 25, 0, 25])
    plt.grid(True)
    return plt
plt = runplt()
X = [[6], [8], [10], [14], [18]]
y = [[7], [9], [13], [17.5], [18]]
plt.plot(X, y, 'k.')
plt.show()

上图中，’x’轴表示匹萨直径，’y’轴表示匹萨价格。能够看出，匹萨价格与其直径正相关，这与我们的
日常经验也比较吻合，自然是越大越贵。下面我们就用scikit-learn来构建模型。

from sklearn import linear_model        #表示，可以调用sklearn中的linear_model模块进行线性回归。
import numpy as np
model = linear_model.LinearRegression()
model.fit(X, y)
display(model.intercept_)  #截距
display(model.coef_)  #线性模型的系数
a = model.predict([[12]])
# a[0][0]
print("预测一张12英寸匹萨价格：{:.2f}".format(model.predict([[12]])[0][0]))

array([ 1.96551724])



array([[ 0.9762931]])


预测一张12英寸匹萨价格：13.68

调用sklearn.linear_model.LinearRegression()所需参数：
* fit_intercept : 布尔型参数，表示是否计算该模型截距。可选参数。
* normalize : 布尔型参数，若为True，则X在回归前进行归一化。可选参数。默认值为False。
* copy_X : 布尔型参数，若为True，则X将被复制；否则将被覆盖。 可选参数。默认值为True。
* n_jobs : 整型参数，表示用于计算的作业数量；若为-1，则用所有的CPU。可选参数。默认值为1。

线性回归fit函数用于拟合输入输出数据，调用形式为model.fit(X,y, sample_weight=None)：
• X : X为训练向量；
• y : y为相对于X的目标向量；
• sample_weight : 分配给各个样本的权重数组，一般不需要使用，可省略。
注意：X，y 以及model.fit()返回的值都是2-D数组，如：a= [ [ 0] ]

一元线性回归假设解释变量和响应变量之间存在线性关系；这个线性模型所构成的空间是一个超平面（hyperplane）。超平面是n维欧氏空间中余维度等于一的线性子空间，**如平面中的直线、空间中的
平面等，总比包含它的空间少一维。在一元线性回归中，一个维度是响应变量，另一个维度是解释变
量，总共两维。因此，其超平面只有一维，就是一条线。**

上述代码中sklearn.linear_model.LinearRegression类是一个估计器（estimator）。**估计
器依据观测值来预测结果**。在scikit-learn里面，所有的估计器都带有**fit()和predict()方
法。fit()用来分析模型参数，predict()是通过**fit()算出的模型参数构成的模型，对解释变量
进行预测获得的值。因为所有的估计器都有这两种方法，所有scikit-learn很容易实验不同的模型。

LinearRegression类的fit()方法学习下面的一元线性回归模型：

$$y = α + βx$$

y表示响应变量的预测值，本例指匹萨价格预测值， x是解释变量，本例指匹萨直径。截距$\alpha$和系数$\beta$ 是线性回归模型最关心的事情。下图中的直线就是匹萨直径与价格的线性关系。用这个模型，你可以计算不同直径的价格，8英寸$7.33，20英寸$18.75。

plt = runplt()
plt.plot(X, y, 'k.')
X2 = [[0], [10], [14], [25]]
model = linear_model.LinearRegression()
model.fit(X,y)
y2 = model.predict(X2)
plt.plot(X2, y2, 'g-')
plt.show()

一元线性回归拟合模型的参数估计常用方法是普通最小二乘法（ordinary least squares ）或线性最小二乘法（linear least squares）。首先，我们定义出拟合成本函数，然后对参数进行数理统计。
附：在线性回归的数学原理blog中将详细细述最小二乘法的前世今生，这里又挖了一坑?，哎，慢慢填吧！

3.2. 带成本函数的模型拟合评估

下图是由若干参数生成的回归直线。如何判断哪一条直线才是最佳拟合呢？

plt = runplt(size=(10,10))
plt.plot(X, y, 'k.')
y3 = [14.25, 14.25, 14.25, 14.25]
y4 = y2*0.5 + 5
model.fit(X[1:-1], y[1:-1])
y5 = model.predict(X2)
plt.plot(X, y, 'k.', label="X, y") 
plt.plot(X2, y2, 'g-.', label="X2 y2")
plt.plot(X2, y3, 'r-.',label="X2, y3")
plt.plot(X2, y4, 'y-.',label="X2, y4")
plt.plot(X2, y5, 'o-', label="X2, y5")
plt.legend()
plt.show()

成本函数（cost function）也叫损失函数（loss function），用来定义模型与观测值的误差。模型预
测的价格与训练集数据的差异称为残差（residuals）或训练误差（training errors）。后面我们会用
模型计算测试集，那时模型预测的价格与测试集数据的差异称为预测误差（prediction errors）或测试误差（test errors）。

附：李航老师的统计学习方法中：将训练误差称为近似误差，将预测误差称为估计误差
* 近似误差：可以理解为对现有训练集的训练误差。
* 估计误差：可以理解为对测试集的测试误差。

模型的残差是训练样本点与线性回归模型的纵向距离，如下图所示

plt = runplt()
plt.plot(X, y, 'k.')
X2 = [[0], [10], [14], [25]]
model = linear_model.LinearRegression()
model.fit(X, y)
y2 = model.predict(X2)
plt.plot(X, y, 'k.')
plt.plot(X2, y2, 'g-')

# 残差预测值
yr = model.predict(X)
#  enumerate 函数可以把一个 list 变成索引-元素对
for idx, x in enumerate(X):
    plt.plot([x, x], [y[idx], yr[idx]], 'r-')
plt.show()

我们可以通过残差之和最小化实现最佳拟合，也就是说模型预测的值与训练集的数据最接近就是**最佳
拟合。对模型的拟合度进行评估的函数称为残差平方和（residual sum of squares）**成本函数。就是
让所有训练数据与模型的残差的平方之和最小化，如下所示：

$$SS_{res} = \sum_{i=1}^n{(y_i - f(x_i))^2}$$

其中，$y_i$ 是观测值，$f(x_i)$ 是预测值。

残差平方和计算如下：

import numpy as np
print('残差平方和:{:.2f}'.format(np.mean((model.predict(X) - y) ** 2)))

残差平方和:1.75

有了成本函数，就要使其最小化获得参数。

3. 3. 解一元线性回归的最小二乘法

通过成本函数最小化获得参数，我们先求系数$\beta$。按照频率论的观点，我们首先需要计算 x的方
差和 x与 y的协方差。

$$var(x) = \frac {\sum_{i=1}^n{(x_i - \overline x)^2 }} {n-1}$$

其中， $\overline x$是直径x 的均值， $x_i$的训练集的第 i个直径样本， n是样本数量。计算如下：

xbar = (6 + 8 + 10 + 14 + 18) / 5
variance = ((6 - xbar)**2 + (8 - xbar)**2 + (10 - xbar)**2 + (14 - xbar)**2 + (18 - xbar)**2) / 4
print(variance)

23.2

Numpy里面有var方法可以直接计算方差，ddof参数是：Delta Degrees of Freedom“：计算中使用的除数是”N-ddof“，其中”N“代表元素的数量。默认情况下，”ddof“为零。设置为1，可得样本方差无偏估计量

import numpy as np
print(np.var([6, 8, 10, 14, 18], ddof=1))

23.2

协方差表示两个变量的总体的变化趋势。如果两个变量的变化趋势一致，也就是说如果其中一个大于
自身的期望值，另外一个也大于自身的期望值，那么两个变量之间的协方差就是正值。 如果两个变
量的变化趋势相反，即其中一个大于自身的期望值，另外一个却小于自身的期望值，那么两个变量之
间的协方差就是负值。如果两个变量不相关，则协方差为0，变量线性无关不表示一定没有其他相关
性。协方差公式如下：

$$cov(x, y) = \frac {\sum_{i=1}^n{(x_i - \overline x)(y_i - \overline y) }} {n-1}$$

其实协方差已经可以表示两个向量之间的关系了，但是会受到向量长度的影响，比如：

虽然两个向量的夹角相等，但是算出来的协方差，除了符号相同外，数值却相差较大，为了解决这个问题，我们把协方差归一化，也就是相关系数。

ybar = (7 + 9 + 13 + 17.5 + 18) / 5
cov = ((6 - xbar) * (7 - ybar) + (8 - xbar) * (9 - ybar) + (10 - xbar) *(13 - ybar) +(14 - xbar) * (17.5 - ybar) + (18 - xbar) * (18 - ybar)) / 4
print(cov)

22.65

Numpy里面有cov方法可以直接计算协方差

import numpy as np
print(np.cov([6, 8, 10, 14, 18], [7, 9, 13, 17.5, 18])[0][1])

22.65

现在有了方差和协方差，就可以计算系数 $\beta$了

算出$\beta$ 后，我们就可以计算$\alpha$ 了:

$$\alpha = y - \beta \overline x$$

将前面的数据带入公式就可以求出α 了:

$$α = 12.9 − 0.9762931034482758 × 11.2 = 1.9655172413793114$$

这样就通过最小化成本函数求出模型参数了。把匹萨直径带入方程就可以求出对应的价格了，如11英
寸直径价格$12.70，18英寸直径价格$19.54。

3.4 模型评估

前面我们用学习算法对训练集进行估计，得出了模型的参数。如何评价模型在现实中的表现呢？**现在
让我们假设有另一组数据，作为测试集进行评估。**

有些度量方法可以用来评估预测效果，我们用R方（r-squared）评估匹萨价格预测的效果。R方也叫
确定系数（coefficient of determination），表示模型对现实数据拟合的程度。计算R方的方法有几
种。一元线性回归中R方等于皮尔逊积矩相关系数（Pearson product moment correlation coefficient
或Pearson’s r）的平方

这种方法计算的R方一定介于0～1之间的正数。其他计算方法，包括scikit-learn中的方法，不是用皮
尔逊积矩相关系数的平方计算的，因此当模型拟合效果很差的时候R方会是负值。下面我们用scikitlearn方法来计算R方。

**附:**R方这里存疑，我不太了解，不确定R方（r-squared）可以评估匹萨价格预测的效果

首先，我们需要计算样本总体平方和， $\overline y$是价格 的均值， $y_i$的训练集的第i 个价格样本，n 是样本数
量。

R方是0.6620说明测试集里面过半数的价格都可以通过模型解释。现在，让我们用scikit-learn来验证
一下。LinearRegression的score方法可以计算R方：

# 测试集
X_test = [[8], [9], [11], [16], [12]]
y_test = [[11], [8.5], [15], [18], [11]]
model = linear_model.LinearRegression()
model.fit(X, y)
model.score(X_test, y_test)

0.6620052929422553