机器学习中的回归分析

通过回归分析预测连续变量

from IPython.display import Image
%matplotlib inline

1.介绍线性回归

线性回归的目标是：对一个或者多个特征与连续型的目标变量之间的关系进行建模。

1.1简单线性回归

对于单个的特征，单变量线性回归表示如下：
$$y=w_0+w_{1}x$$

其中，$w_0$代表的是$y$轴截距，$w_1$代表特征的权重系数。目标就是：学习线性方程的权重，从而描述特征变量与目标变量之间的关系，然后利用这些权重来预测不同训练集的新的特征变量所对应的目标变量。

# 线性回归几何表示
Image(filename='images/10_01.png', width=500) 

最佳拟合直线也叫作回归线，预测结果与真实结果之间的垂直线条叫做偏移或者残差。

1.2多变量线性回归

多变量线性回归表示如下：
$$y=w_0{x}_0+w_1{x}_1+\dots +w_m{x}_m=\sum _{i=0}^n{w}_i{x}_i=w^{T}x$$

这里的$w_0$为$y$轴截距，同时$x_0$为1.

Image(filename='images/10_15.png', width=500) 

2.探索房屋数据

2.1将数据加载成一个DataFrame

数据下载地址： https://archive.ics.uci.edu/ml/datasets/Housing

数据集合有506个样本，特征描述:

1. CRIM      各个城镇人均犯罪率；
2. ZN       宅基地面积超过25000平方米的占地比例；
3. INDUS     城镇非零售经营面积比例；
4. CHAS      如果区域边界为河流，取值为1，否则取值为0；
5. NOX       一氧化碳浓度；
6. RM        每一户平均房间数；
7. AGE       1940年之前建成的业主自住单位的比例；
8. DIS       到波士顿五个就业中心的加权距离；
9. RAD       放射状公路通达性指数；
10. TAX      每万元的全价值物业税税率；
11. PTRATIO   按镇划分的师生比例；
12. B       1000(BK-0.63)**2，其中，BK代表的是城镇中非裔美国人的比例；
13. LSTAT     地位较低的人口百分比；
14. MEDV     自住房屋价值中位数；

# 载入数据
import pandas as pd

df = pd.read_csv('./housing.data.txt',
                 header=None,
                 sep='\s+')   # \s+ 意思就是至少有一个空白字符存在

df.columns = ['CRIM', 'ZN', 'INDUS', 'CHAS', 
              'NOX', 'RM', 'AGE', 'DIS', 'RAD', 
              'TAX', 'PTRATIO', 'B', 'LSTAT', 'MEDV']
df.head()

	CRIM	ZN	INDUS	CHAS	NOX	RM	AGE	DIS	RAD	TAX	PTRATIO	B	LSTAT	MEDV
0	0.00632	18.0	2.31	0	0.538	6.575	65.2	4.0900	1	296.0	15.3	396.90	4.98	24.0
1	0.02731	0.0	7.07	0	0.469	6.421	78.9	4.9671	2	242.0	17.8	396.90	9.14	21.6
2	0.02729	0.0	7.07	0	0.469	7.185	61.1	4.9671	2	242.0	17.8	392.83	4.03	34.7
3	0.03237	0.0	2.18	0	0.458	6.998	45.8	6.0622	3	222.0	18.7	394.63	2.94	33.4
4	0.06905	0.0	2.18	0	0.458	7.147	54.2	6.0622	3	222.0	18.7	396.90	5.33	36.2

2.2通过可视化来探索数据集中的一些重要特性

import matplotlib.pyplot as plt
# 这里从mlxtend工具包中导入散点图矩阵工具
from mlxtend.plotting import scatterplotmatrix

这里通过scatterplot matrix来可视化成对特征之间的相关性。这里由于页面尺寸的限制，仅选择了5个特征

cols = ['LSTAT', 'INDUS', 'NOX', 'RM', 'MEDV']

scatterplotmatrix(df[cols].values, figsize=(10, 8), 
                  names=cols, alpha=0.5)
plt.tight_layout()
#plt.savefig('images/10_03.png', dpi=300)
plt.show()

通过散点图矩阵，可以很方便地看出数据分布情况以及其是否包含异常值，可以看出RM和和MEDV之间呈现出线性关系，同时MEDV看起来呈现出正态分布形式，同时包含一定的异常值。

cols = ['CRIM', 'ZN', 'CHAS', 
               'AGE', 'MEDV']

scatterplotmatrix(df[cols].values, figsize=(10, 8), 
                  names=cols, alpha=0.5)
plt.tight_layout()
#plt.savefig('images/10_03.png', dpi=300)
plt.show()

训练一个线性回归模型并未要求特征变量和目标变量满足正态分布。正态分布仅为某些数据统计和假设检验的要求。

2.3使用相关性矩阵来反映特征之间的相关性

皮尔森矩积相关系数： 度量了成对特征之间的线性依赖关系。取值在正负一之间。对特征$x$和$y$：
$$r=\frac{{\sum }_{i=1}^n\left[\left(x^{(i)}-{\mu }_x\right)\left(y^{(i)}-{\mu }_y\right)\right]}{\sqrt{{\sum }_{i=1}^n{\left(x^{(i)}-{\mu }_x\right)}^2}\sqrt{{\sum }_{i=1}^n{\left(y^{(i)}-{\mu }_y\right)}^2}}=\frac{{\sigma }_{xy}}{{\sigma }_x{\sigma }_y}$$

这里的$\mu$代表的是均值，${\sigma }_{xy}$代表的是特征协方差，${\sigma }_x$和${\sigma }_y$代表的是特征的标准差

**标准化特征的协方差和相关性：**可以证明，一对标准化特征之间的协方差等于他们之间的线性相关系数。

这里首先对特征$x$和$y$进行特征标准化处理：
$$x^{\prime }=\frac{x-{\mu }_x}{{\sigma }_x},y^{\prime }=\frac{y-{\mu }_y}{{\sigma }_y}$$

计算协方差：

$${\sigma }_{xy}=\frac{1}{n}\sum _i^n\left(x^{(i)}-{\mu }_x\right)\left(y^{(i)}-{\mu }_y\right)$$

标准化之后，特征的均值为零，标准差为1.因此：

$${\sigma }_{xy}^{\prime }=\frac{1}{n}\sum _i^n\left(x^{\prime (i)}-0\right)\left(y^{\prime (i)}-0\right)$$

将上式中标准化之后的特征用标准化后的形式替换：

$$\begin{array}{rl}{\sigma }_{xy}^{\prime }& =\frac{1}{n}\sum _i^n\left(\frac{x-{\mu }_x}{{\sigma }_x}\right)\left(\frac{y-{\mu }_y}{{\sigma }_y}\right)\\ {\sigma }_{xy}^{\prime }& =\frac{1}{n\cdot {\sigma }_x{\sigma }_y}\sum _i^n\left(x^{(i)}-{\mu }_x\right)\left(y^{(i)}-{\mu }_y\right)\end{array}$$

因此得到：

$${\sigma }_{xy}^{\prime }=\frac{{\sigma }_{xy}}{{\sigma }_x{\sigma }_y}$$

下面将使用Numpy的corrcoef函数在上述散点图矩阵的5个特征基础上，结合Mlxtend的heatmap函数来绘制热力图。

cols = ['CRIM', 'ZN', 'INDUS', 'CHAS', 
              'NOX', 'RM', 'AGE', 'DIS', 'RAD', 
              'TAX', 'PTRATIO', 'B', 'LSTAT', 'MEDV']

import numpy as np
from mlxtend.plotting import heatmap


cm = np.corrcoef(df[cols].values.T)
hm = heatmap(cm, row_names=cols, column_names=cols, figsize=(15,15))

# plt.savefig('images/10_04.png', dpi=300)
plt.show()

为了拟合线性回归模型，我们感兴趣的是那些与目标变量MEDV有高相关性的特征。从heatmap中可以看到，目标变量MEDV与LSTAT变量的相关性

最大(-0.74);但是之前在检查散点图矩阵scatterplotmatrix时，LSTAT和MEDV之间存在明显的非线性关系。另一方面，RM与MEDV的相关性也较高

(0.70)。鉴于我们在散点图中观察到的这两个变量之间的线性关系，RM似乎是一个很好的选择。

3.实现普通最小二乘线性回归模型

使用**Ordinary least squares （OLS）**最小二乘法，有时候也称为线性最小二乘，来估计线性回归模型在具有最小误差平方和的时候的模型参数。

3.1使用梯度下降法求解回归参数

在自适应线性神经元中，其使用的是一个线性激活函数。定义的损失函数为$J(W)$，同时结合梯度下降或者随机梯度下降来实现优化。

Adaline中的损失函数为SSE(残差平方和)，与最小二乘法OLS中使用的损失函数相同：

$$J(w)=\frac{1}{2}\sum _{i=1}^n{\left(y^{(i)}-{\hat{y}}^{(i)}\right)}^2$$

其中，$y^{(i)}$为预测值$W^{T}X$，OLS回归可以视为AdaLine在没有单位阶跃函数的情况下，产生的连续值输出。下面是GradientDescent的实现：

class LinearRegressionGD(object):

    def __init__(self, eta=0.001, n_iter=20):
        self.eta = eta
        self.n_iter = n_iter

    def fit(self, X, y):
        self.w_ = np.zeros(1 + X.shape[1])
        self.cost_ = []

        for i in range(self.n_iter):
            output = self.net_input(X)
            errors = (y - output)
            self.w_[1:] += self.eta * X.T.dot(errors)
            self.w_[0] += self.eta * errors.sum()
            cost = (errors**2).sum() / 2.0
            self.cost_.append(cost)
        return self

    def net_input(self, X):
        return np.dot(X, self.w_[1:]) + self.w_[0]

    def predict(self, X):
        return self.net_input(X)

详细的权重更新过程

# 下面以RM和MEDV为示例，展示回归变量在回归模型中的作用，
X = df[['RM']].values
y = df['MEDV'].values

from sklearn.preprocessing import StandardScaler
# 这里为了使得GD算法更快地收敛，对特征进行标准化

sc_x = StandardScaler()
sc_y = StandardScaler()
X_std = sc_x.fit_transform(X)
y_std = sc_y.fit_transform(y[:, np.newaxis]).flatten()

"""
sklearn中大部分的transformers都希望数据存储在二维数组中，这里使用np.newaxis新增了一个维度。

然后再StandardScaler返回缩放之后的特征后，使用flatten将其转换为一维数组表示形式。
"""

'\nsklearn中大部分的transformers都希望数据存储在二维数组中，这里使用np.newaxis新增了一个维度。\n\n然后再StandardScaler返回缩放之后的特征后，使用flatten将其转换为一维数组表示形式。\n'

lr = LinearRegressionGD()
lr.fit(X_std, y_std)

<__main__.LinearRegressionGD at 0x258c7072288>

# 将损失函数绘制为epochs的函数(基于完整的iterations)，有助于更明确地反映算法收敛情况
plt.plot(range(1, lr.n_iter+1), lr.cost_)
plt.ylabel('SSE')
plt.xlabel('Epoch')
#plt.tight_layout()
#plt.savefig('images/10_05.png', dpi=300)
plt.show()

可以看到GD算法在第5个epoch的时候，收敛了

# 绘制回归直线，查看拟合情况
def lin_regplot(X, y, model):
    plt.scatter(X, y, c='steelblue', edgecolor='white', s=70)
    plt.plot(X, model.predict(X), color='black', lw=2)    
    return 

lin_regplot(X_std, y_std, lr)
plt.xlabel('Average number of rooms [RM] (standardized)')
plt.ylabel('Price in $1000s [MEDV] (standardized)')

#plt.savefig('images/10_06.png', dpi=300)
plt.show()

回归直线大致反映出房价随着房屋数量的增加而呈现出上涨趋势。但很多取值情况下，房间数量并不能很好地解释房价。同时可以看出，若干数据点在y=3

的时候，以此排开，这说明房价可能被削减抑制了。

print('Slope: %.3f' % lr.w_[1]) # 回归直线斜率
print('Intercept: %.3f' % lr.w_[0]) # 回归直线截距

Slope: 0.695
Intercept: -0.000

注意：如果使用的是标准化特征，则从技术上讲，无需更新截距intercept权重，因为这种情况下，y轴截距始终为0。上面是打印结果；

为了将预测结果以原始尺度进行显示，可以直接通过StandardScaler的inverse_transform方法逆向映射

num_rooms_std = sc_x.transform(np.array([[5.0]]))
price_std = lr.predict(num_rooms_std)
print("Price in $1000s: %.3f" % sc_y.inverse_transform(price_std))

Price in $1000s: 10.840

3.2通过sklearn估计回归模型的系数

sklearn中的线性回归也适用于非标准化特征(效果更好)，因为它没有使用基于GD的优化，因此可以跳过标准化步骤：

from sklearn.linear_model import LinearRegression

slr = LinearRegression()
slr.fit(X, y)

y_pred = slr.predict(X)
print('Slope: %.3f' % slr.coef_[0])
print('Intercept: %.3f' % slr.intercept_)        # 这里没有进行标准化，产生的模型系数也不同，但是与前面的拟合图相比，这里是一样的

Slope: 9.102
Intercept: -34.671

lin_regplot(X, y, slr)
plt.xlabel('Average number of rooms [RM]')
plt.ylabel('Price in $1000s [MEDV]')

#plt.savefig('images/10_07.png', dpi=300)
plt.show()

线性回归的解析解：

$$\mathbf{w}={\left({\mathbf{X}}^T\mathbf{X}\right)}^{-1}{\mathbf{X}}^T\mathbf{y}$$

$$\begin{array}{c}{\nabla }_{\mathbf{w}}{\mathrm{M}\mathrm{S}\mathrm{E}}_{\text{train }}=0\\ \Rightarrow {\nabla }_{\mathbf{w}}\frac{1}{m}{\Vert {\hat{\mathbf{y}}}^{(\text{train })}-{\mathbf{y}}^{(\text{train })}\Vert }_2^2=0\\ \Rightarrow \frac{1}{m}{\nabla }_{\mathbf{w}}{\Vert {\mathbf{X}}^{(\text{train })}\mathbf{w}-{\mathbf{y}}^{(\text{train })}\Vert }_2^2=0\end{array}$$

$$\begin{array}{c}\Rightarrow {\nabla }_{\mathbf{w}}{\left({\mathbf{X}}^{(\text{train })}\mathbf{w}-{\mathbf{y}}^{(\text{train })}\right)}^{\top }\left({\mathbf{X}}^{(\text{train })}\mathbf{w}-{\mathbf{y}}^{(\text{train })}\right)=0\\ \Rightarrow {\nabla }_{\mathbf{w}}\left({\mathbf{w}}^{\top }{\mathbf{X}}^{(\mathrm{train})\top }{\mathbf{X}}^{(\mathrm{train})}\mathbf{w}-2{\mathbf{w}}^{\top }{\mathbf{X}}^{(\mathrm{train})\top }{\mathbf{y}}^{(\text{train })}+{\mathbf{y}}^{(\mathrm{train})\top }{\mathbf{y}}^{(\text{train })}\right)=0\\ \Rightarrow 2{\mathbf{X}}^{(\text{train })\top }{\mathbf{X}}^{(\text{train })}\mathbf{w}-2{\mathbf{X}}^{(\mathrm{train})\top }{\mathbf{y}}^{(\text{train })}=0\\ \Rightarrow \mathbf{w}={\left({\mathbf{X}}^{(\mathrm{train})\top }{\mathbf{X}}^{(\mathrm{train})}\right)}^{-1}{\mathbf{X}}^{(\mathrm{train})\top }{\mathbf{y}}^{(\text{train })}\end{array}$$

# python实现
# adding a column vector of "ones"
Xb = np.hstack((np.ones((X.shape[0], 1)), X))
w = np.zeros(X.shape[1])

z = np.linalg.inv(np.dot(Xb.T, Xb))
w = np.dot(z, np.dot(Xb.T, y))

print('Slope: %.3f' % w[1])
print('Intercept: %.3f' % w[0])

Slope: 9.102
Intercept: -34.671

上述方法的优点：保证找到最优解。但当数据集非常大，则上式求解逆矩阵的代价就相当高。

此外，如果包含训练样本的矩阵是奇异矩阵，不可逆，则也会面临求解困难。因此有时候会更偏向于迭代求解方法。

4.使用随机样本一致性RANSAC训练更加鲁棒的回归模型

线性回归模型会受到异常值的严重影响，而且在一些情况下，数据中非常小的子集也可能对模型系数产生巨大影响。可以通过统计检验的方法来检测异常值。

这里采用随机样本一致性RANSAC算法来作为另一种更具鲁棒性的回归方法。

RANSAC算法步骤总结：

1.选择随机数量的样本，作为内点inliers并训练模型；

2.根据拟合模型测试其他所有的数据点，并将在用户定义的公差tolerance范围内的点都添加到内点inliers中；

3.使用所有的inliers重新训练模型；

4.评估模型相对于内点inliers的误差；

5.如果模型性能达到了用户定义的某个阈值或者达到了最大迭代次数，则终止算法。否则返回步骤1；

# RANSAC在sklearn中已经做了实现
from sklearn.linear_model import RANSACRegressor

ransac = RANSACRegressor(LinearRegression(), 
                         max_trials=100,   # 最大迭代次数
                         min_samples=50,   # 随机选择的最少样本数
                         loss='absolute_loss',  # 使用绝对误差
                         residual_threshold=5.0,  # 残差阈值设定为5，这里当预测结果与真实值之间的垂直偏移在5以内，则将当前样本视为inliers
                         random_state=0)


ransac.fit(X, y)

"""
拟合了RANSAC之后，就可以得到内点和异常值，这里通过绘制拟合直线展示
"""

inlier_mask = ransac.inlier_mask_
outlier_mask = np.logical_not(inlier_mask)

line_X = np.arange(3, 10, 1)
line_y_ransac = ransac.predict(line_X[:, np.newaxis])
plt.scatter(X[inlier_mask], y[inlier_mask],
            c='steelblue', edgecolor='white', 
            marker='o', label='Inliers')
plt.scatter(X[outlier_mask], y[outlier_mask],
            c='limegreen', edgecolor='white', 
            marker='s', label='Outliers')
plt.plot(line_X, line_y_ransac, color='black', lw=2)   
plt.xlabel('Average number of rooms [RM]')
plt.ylabel('Price in $1000s [MEDV]')
plt.legend(loc='upper left')

#plt.savefig('images/10_08.png', dpi=300)
plt.show()

默认情况下，sklearn通过MAD来选择上述的“inlier threshold”。其中，MAD代表的是目标值y的中值绝对偏差median absolute deviation。

但实际上，“inlier threshold”的选择需要针对具体问题，这是RANSAC的一个缺陷。相关论文指出可以自动选择内点阈值。

# 拟合结果与前面的有一定区别，但使用RANSAC降低了数据集中异常值的潜在影响
print('Slope: %.3f' % ransac.estimator_.coef_[0])
print('Intercept: %.3f' % ransac.estimator_.intercept_)

Slope: 10.735
Intercept: -44.089

5.评价线性回归模型的性能

from sklearn.model_selection import train_test_split

X = df.iloc[:, :-1].values
y = df['MEDV'].values

# 训练集，测试集划分
X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.3, random_state=0)

slr = LinearRegression()

slr.fit(X_train, y_train)

y_train_pred = slr.predict(X_train)
y_test_pred = slr.predict(X_test)

import numpy as np
import scipy as sp

ary = np.array(range(100000))

%timeit np.linalg.norm(ary)

138 µs ± 7.65 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

%timeit sp.linalg.norm(ary)

141 µs ± 4.19 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

%timeit np.sqrt(np.sum(ary**2))

167 µs ± 1.19 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

模型使用了多个解释变量，所以不能将线性回归线(超平面)可视化在二维图中，这里选择绘制残差(实际值和预测值之间的差异或垂直距离)与预测值之间的关

系，以评估回归模型。残差图是诊断回归模型的常用图形工具。可以帮助检测非线性和离群值，并检查误差是否随机分布。使用以下代码绘制残差图：

plt.scatter(y_train_pred,  y_train_pred - y_train,
            c='steelblue', marker='o', edgecolor='white',
            label='Training data')
plt.scatter(y_test_pred,  y_test_pred - y_test,
            c='limegreen', marker='s', edgecolor='white',
            label='Test data')

plt.xlabel('Predicted values')
plt.ylabel('Residuals')
plt.legend(loc='upper left')
plt.hlines(y=0, xmin=-10, xmax=50, color='black', lw=2)
plt.xlim([-10, 50])
plt.tight_layout()

# plt.savefig('images/10_09.png', dpi=300)
plt.show()

在完全准确预测的情况下，残差正好为零，但是这种情况在实际应用中基本遇不到。

对于优秀的回归模型，其误差应该是随机分布的，也就是随机分布在中心线周围。如果残差图呈现出固定的规律模式，则说明模型无法捕获到一些已经

泄漏到残差中的解释性信息。

同时，残差图也可以用于检测离群点，离群点由偏离中心线较大的点表示。

模型性能的另一个有用的指标为MSE，它是我们试图最小化的误差平方和的平均值。

MSE在用于对比不同的回归模型或者通过网络搜索和交叉验证调参的时候非常有用，因为它根据样本大小归一化了SSE：

$$MSE=\frac{1}{n}\sum _{i=1}^n{\left(y^{(i)}-{\hat{y}}^{(i)}\right)}^2$$

下面计算了模型在测试集和训练集上的MSE：

from sklearn.metrics import r2_score
from sklearn.metrics import mean_squared_error

print('MSE train: %.3f, test: %.3f' % (
        mean_squared_error(y_train, y_train_pred),
        mean_squared_error(y_test, y_test_pred)))

"""
可以看到，模型在训练集上的MSE为19.9958，而测试集上的MSE为27.196，这说明模型有一定的过拟合。

但是MSE与模型分类精度相比，MSE是无界的。或者说对于MSE的解释依赖于数据集合特征缩放。

例如：将房价表示为1000的倍数(使用K表示后缀)，与未使用缩放特征的模型相比，相同的模型可能会产生较低的MSE。($10K-15K)^2 < ($10000-$15000)^2
"""
print('R^2 train: %.3f, test: %.3f' % (
        r2_score(y_train, y_train_pred),
        r2_score(y_test, y_test_pred)))

MSE train: 20.926, test: 28.876
R^2 train: 0.753, test: 0.653

可以看到，模型在训练集上的MSE为19.9958，而测试集上的MSE为27.196，这说明模型有一定的过拟合。

但是MSE与模型分类精度相比，MSE是无界的。或者说对于MSE的解释依赖于数据集合特征缩放。

例如：将房价表示为1000的倍数(使用K表示后缀)，与未使用缩放特征的模型相比，相同的模型可能会产生较低的MSE。$($10K-$15K{)}^2$ < $($10000-$15000{)}^2$

因此,有时候决定系数 cofficient of determination 可能更加有用,它可以理解为MSE的标准化版本，以便更好地解释模型性能。

$R^2$定义如下：

$$R^2=1-\frac{SSE}{SST}$$

其中，SSE为误差平方和，SST为目标变量的方差

$$SST=\sum _{i=1}^n{\left(y^{(i)}-{\mu }_y\right)}^2$$

$R^2$实际上是MSE的一个缩放版本：

$$\begin{array}{rl}{R}^2& =1-\frac{SSE}{SST}\\ & =\frac{\frac{1}{n}{\sum }_{i=1}^n{\left(y^{(i)}-{\hat{y}}^{(i)}\right)}^2}{\frac{1}{n}{\sum }_{i=1}^n{\left(y^{(i)}-{\mu }_y\right)}^2}\\ & =1-\frac{MSE}{\mathrm{Var}(y)}\end{array}$$

对于训练集，$R^2$被限定在0和1之间，但是对于测试集，它有可能取到负值。

如果$R^2=1$,说明模型对数据的拟合非常完美，同时等价于$MSE=0$。

这里对模型的评估上，测试集上的$R^2=0.765$,在测试集上$R^2=0.673$。确实有点过拟合了

6.使用正则化方法进行回归分析

正则化是通过添加对模型权重的惩罚来抑制模型复杂度，从而环节过拟合的一种方法。

最流行的正则化回归方法有：岭回归、最小绝对收缩least absolute shrinkage和套索回归，以及弹性网络elastic Net。

对于岭回归：是在误差平方和之后加上模型权重的L2正则项：

$$J(\mathbf{w}{)}_{\text{Ridge }}=\sum _{i=1}^n{\left(y^{(i)}-{\hat{y}}^{(i)}\right)}^2+\lambda \Vert w{\Vert }_2^2$$

其中，对于正则项：

$$L2:\lambda \Vert w{\Vert }_2^2=\lambda \sum _{j=1}^m{w}_j^2$$

通过增大超参数$\lambda$的值，增加了对模型权重$w$的惩罚程度。

注意，正则化项中不包含$y$轴截距$W_0$。

另一种常用的方法是：套索回归，归因于L1正则项，很多的权重值都变成了0，因此套索回归也经常被用于特征选择。

L1正则化特征选择，程序示例。

$$J(w{)}_{LASSO}=\sum _{i=1}^n{\left(y^{(i)}-{\hat{y}}^{(i)}\right)}^2+\lambda \Vert w{\Vert }_1$$

其中，L1惩罚定义为模型权重的绝对值之和：

$$L1:\lambda \Vert w{\Vert }_1=\lambda \sum _{j=1}^m\left|w_j\right|$$

Lasso算法的局限性：

当$m>n$的时候，它最多只能选择$n$个特征，其中$n$是训练样本的数量。这在特征选择的某些应用中可能是不希望的。然而，在实践中，套索的这一特性通常是一

种优势，因为它避免了饱和模型。如果训练样本的数量等于特征的数量，则会出现模型饱和，这是过度参数化的一种形式。因此，饱和模型总是可以很好地

拟合训练数据，但只是一种插值形式，因此不能很好地推广。

L1和L2惩罚的一种折中形式：弹性网络

其中，L1惩罚用于生成稀疏解，L2惩罚用于选择大于样本个数的$n$个特征，(当出现m>n)的情况：

$$J(\mathbf{w}{)}_{\text{ElasticNet }}=\sum _{i=1}^n{\left(y^{(i)}-{\hat{y}}^{(i)}\right)}^2+{\lambda }_1\sum _{j=1}^m{w}_j^2+{\lambda }_2\sum _{j=1}^m\left|w_j\right|$$

其中，对于正则项系数$\lambda$的选择，通常可以结果K-fold交叉验证来实现择优。

# sklearn中Lasso的实现
from sklearn.linear_model import Lasso

lasso = Lasso(alpha=0.1)
lasso.fit(X_train, y_train)

y_train_pred = lasso.predict(X_train)
y_test_pred = lasso.predict(X_test)
print(lasso.coef_)

[-0.11311792  0.04725111 -0.03992527  0.96478874 -0.          3.72289616
 -0.02143106 -1.23370405  0.20469    -0.0129439  -0.85269025  0.00795847
 -0.52392362]

print('MSE train: %.3f, test: %.3f' % (
        mean_squared_error(y_train, y_train_pred),
        mean_squared_error(y_test, y_test_pred)))
print('R^2 train: %.3f, test: %.3f' % (
        r2_score(y_train, y_train_pred),
        r2_score(y_test, y_test_pred)))

MSE train: 20.926, test: 28.876
R^2 train: 0.753, test: 0.653

sklearn中的Ridge regression实现：

from sklearn.linear_model import Ridge
ridge = Ridge(alpha=1.0)

ridge.fit(X_train, y_train)

y_train_pred = ridge.predict(X_train)
y_test_pred = ridge.predict(X_test)
print(ridge.coef_)

print('MSE train: %.3f, test: %.3f' % (
        mean_squared_error(y_train, y_train_pred),
        mean_squared_error(y_test, y_test_pred)))
print('R^2 train: %.3f, test: %.3f' % (
        r2_score(y_train, y_train_pred),
        r2_score(y_test, y_test_pred)))

[-1.18308575e-01  4.61259764e-02 -2.08626416e-02  2.45868617e+00
 -8.25958494e+00  3.89748516e+00 -1.79140171e-02 -1.39737175e+00
  2.18432298e-01 -1.16338128e-02 -9.31711410e-01  7.26996266e-03
 -4.94046539e-01]
MSE train: 20.145, test: 27.762
R^2 train: 0.762, test: 0.667

sklearn中LASSO regression 的实现：

from sklearn.linear_model import Lasso
lasso = Lasso(alpha=1.0)

lasso.fit(X_train, y_train)

y_train_pred = lasso.predict(X_train)
y_test_pred = lasso.predict(X_test)
print(lasso.coef_)

print('MSE train: %.3f, test: %.3f' % (
        mean_squared_error(y_train, y_train_pred),
        mean_squared_error(y_test, y_test_pred)))
print('R^2 train: %.3f, test: %.3f' % (
        r2_score(y_train, y_train_pred),
        r2_score(y_test, y_test_pred)))

[-0.06586193  0.04832933 -0.          0.         -0.          0.86898466
  0.01217999 -0.75109378  0.2000743  -0.01395062 -0.84602363  0.00668818
 -0.73266568]
MSE train: 24.717, test: 32.345
R^2 train: 0.708, test: 0.612

sklearn中Elastic Net regression 的实现：

from sklearn.linear_model import ElasticNet
elanet = ElasticNet(alpha=1.0, l1_ratio=0.5)

elanet.fit(X_train, y_train)

y_train_pred = elanet.predict(X_train)
y_test_pred = elanet.predict(X_test)
print(elanet.coef_)

print('MSE train: %.3f, test: %.3f' % (
        mean_squared_error(y_train, y_train_pred),
        mean_squared_error(y_test, y_test_pred)))
print('R^2 train: %.3f, test: %.3f' % (
        r2_score(y_train, y_train_pred),
        r2_score(y_test, y_test_pred)))

[-0.08344408  0.05179376 -0.01702468  0.         -0.          0.90890973
  0.01218953 -0.83010765  0.23558231 -0.01502425 -0.84881663  0.00687826
 -0.72504946]
MSE train: 24.381, test: 31.874
R^2 train: 0.712, test: 0.617

7.将线性回归转化为曲线多项式回归

前面的简单回归中都是假设特征与目标变量之间为线性关系，对于不满足线性假设的场合，接下来将使用多项式回归模型

形如：

$$y=w_0+w_{1}x+w_2{x}^2+\dots +w_d{x}^d$$

这里的$d$代表的是多项式的阶。

尽管多项式回归可以被用于处理非线性关系，但由于线性回归系数$w$,他任然被认为是多元线性回归模型。

X = np.array([258.0, 270.0, 294.0, 
              320.0, 342.0, 368.0, 
              396.0, 446.0, 480.0, 586.0])\
             [:, np.newaxis]

y = np.array([236.4, 234.4, 252.8, 
              298.6, 314.2, 342.2, 
              360.8, 368.0, 391.2,
              390.8])

from sklearn.preprocessing import PolynomialFeatures

lr = LinearRegression()
pr = LinearRegression()

# 设定为2阶
quadratic = PolynomialFeatures(degree=2)
X_quad = quadratic.fit_transform(X)

# 拟合原始线性形式的特征，训练简单线性回归模型用于对比
lr.fit(X, y)
X_fit = np.arange(250, 600, 10)[:, np.newaxis]
y_lin_fit = lr.predict(X_fit)

# 对变换后的特征拟合多元回归模型进行多项式回归
pr.fit(X_quad, y)
y_quad_fit = pr.predict(quadratic.fit_transform(X_fit))

# 结果可视化
plt.scatter(X, y, label='Training points')
plt.plot(X_fit, y_lin_fit, label='Linear fit', linestyle='--')
plt.plot(X_fit, y_quad_fit, label='Quadratic fit')

plt.xlabel('Explanatory variable')
plt.ylabel('Predicted or known target values')
plt.legend(loc='upper left')

plt.tight_layout()
#plt.savefig('images/10_11.png', dpi=300)
plt.show()

可以看出，多项式回归对于特征和目标变量之间的关系把握更加紧密一些。

y_lin_pred = lr.predict(X)
y_quad_pred = pr.predict(X_quad)

print('Training MSE linear: %.3f, quadratic: %.3f' % (
        mean_squared_error(y, y_lin_pred),
        mean_squared_error(y, y_quad_pred)))
print('Training R^2 linear: %.3f, quadratic: %.3f' % (
        r2_score(y, y_lin_pred),
        r2_score(y, y_quad_pred)))

Training MSE linear: 569.780, quadratic: 61.330
Training R^2 linear: 0.832, quadratic: 0.982

7.1住房数据集中的非线性关系建模

下面将使用二阶和3阶多项式对房屋价格和LSTAT之间的关系进行建模，并将其与线性拟合进行比较。

X = df[['LSTAT']].values
y = df['MEDV'].values

regr = LinearRegression()

# 通过非线性特征变换，构造非线性特征
quadratic = PolynomialFeatures(degree=2)
cubic = PolynomialFeatures(degree=3)

X_quad = quadratic.fit_transform(X)
X_cubic = cubic.fit_transform(X)

# fit features
X_fit = np.arange(X.min(), X.max(), 1)[:, np.newaxis]

regr = regr.fit(X, y)
y_lin_fit = regr.predict(X_fit)
linear_r2 = r2_score(y, regr.predict(X))

regr = regr.fit(X_quad, y)
y_quad_fit = regr.predict(quadratic.fit_transform(X_fit))
quadratic_r2 = r2_score(y, regr.predict(X_quad))

regr = regr.fit(X_cubic, y)
y_cubic_fit = regr.predict(cubic.fit_transform(X_fit))
cubic_r2 = r2_score(y, regr.predict(X_cubic))


# 结果可视化
plt.scatter(X, y, label='Training points', color='lightgray')

plt.plot(X_fit, y_lin_fit, 
         label='Linear (d=1), $R^2=%.2f$' % linear_r2, 
         color='blue', 
         lw=2, 
         linestyle=':')

plt.plot(X_fit, y_quad_fit, 
         label='Quadratic (d=2), $R^2=%.2f$' % quadratic_r2,
         color='red', 
         lw=2,
         linestyle='-')

plt.plot(X_fit, y_cubic_fit, 
         label='Cubic (d=3), $R^2=%.2f$' % cubic_r2,
         color='green', 
         lw=2, 
         linestyle='--')

plt.xlabel('% lower status of the population [LSTAT]')
plt.ylabel('Price in $1000s [MEDV]')
plt.legend(loc='upper right')

#plt.savefig('images/10_12.png', dpi=300)
plt.show()

7.2对数变换

添加更多的多项式特征可能会造成模型的过拟合。因此在实践中，通常会根据单独的测试集数据评估模型的性能，从而估计模型泛化能力。

同时多项式特征并不总是建模非线性关系的最佳选择。

例如：对于MEDV和LSTAT的散点图矩阵，可能会假设LSTAT特征变量和MEDV的平方根的对数变换可以将数据投影到适合线性回归拟合的线性特征空间。

或者，从散点图矩阵中可以大致看出，LSTAT与MEDV之间的关系类似于指数函数：

$$f(x)=e^{-x}$$

由于指数函数的自然对数是一条直线，因此可以进行一下变换：

$$\log (f(x))=-x$$

预处理代码如下：

X = df[['LSTAT']].values
y = df['MEDV'].values

# 特征变换
X_log = np.log(X)
y_sqrt = np.sqrt(y)

# fit features
X_fit = np.arange(X_log.min()-1, X_log.max()+1, 1)[:, np.newaxis]

regr = regr.fit(X_log, y_sqrt)
y_lin_fit = regr.predict(X_fit)
linear_r2 = r2_score(y_sqrt, regr.predict(X_log))

# 结果可视化
plt.scatter(X_log, y_sqrt, label='Training points', color='lightgray')

plt.plot(X_fit, y_lin_fit, 
         label='Linear (d=1), $R^2=%.2f$' % linear_r2, 
         color='blue', 
         lw=2)

plt.xlabel('log(% lower status of the population [LSTAT])')
plt.ylabel('$\sqrt{Price \; in \; \$1000s \; [MEDV]}$')
plt.legend(loc='lower left')

plt.tight_layout()
#plt.savefig('images/10_13.png', dpi=300)
plt.show()

将特征转换到对数空间，对目标变量取平方根之后，对于原本无法使用回归直线捕获的关系可以直接通过一条线性回归线进行捕捉。同时优于上面的多项式回归。

8.使用随机森林处理非线性关系

随机森林是多个决策树的集成，可以理解为分段线性函数的和，而不是全局线性和多项式回归模型。

通过决策树算法，将输入空间细分为更小的区域，从而更易于管理。

8.1决策树回归

决策树算法的一个优点：在处理非线性数据的时候，无需进行任何变换，因为决策树一次性考虑一个特征，而不是特征的加权组合。

同时，决策树算法也无需进行标准化或者归一化处理，因为标准化或者归一化不改变特征的信息增益等指标。

当使用决策树分类的时候，将信息熵定义为信息的不纯度，从而确定特征分裂对于信息增益的最大化，定义如下：

$$IG\left(D_p,x_i\right)=I\left(D_p\right)-\frac{N_{\text{left }}}{N_p}I\left(D_{\text{left }}\right)-\frac{N_{\text{right }}}{N_p}I\left(D_{\text{right }}\right)$$

这里的$x_i$为当前正要分裂的特征，$N_p$是父节点中训练样本的数量，$I$为熵计算函数，$D_p$为父节点中训练样本子集；

$D_{left}$和$D_{right}$为分裂之后，左右子节点中的训练样本子集；

目标是发现特征分裂之后，最大化信息增益的分裂方式。也就是要找到分裂之后，造成信息不纯度减少最多的特征。除开基尼不纯度等指标，这里为了

适用于连续变量，所以结合MSE定义了结点$t$的不纯度指标：

$$I(t)=MSE(t)=\frac{1}{N_t}\sum _{i\in D_t}{\left(y^{(i)}-{\hat{y}}_t\right)}^2$$

其中，$N_t$为结点$t$中的训练样本数；$D_t$为结点$t$的训练子集；$y^{(i)}$代表的是真实目标变量，${\hat{y}}_t$为预测结果；

$${\hat{y}}_t=\frac{1}{N_t}\sum _{i\in D_t}{y}^{(i)}$$

在决策树回归的背景下，MSE通常被称为whthin-node variance，这就是将方差降低作为分裂标准的原因。

from sklearn.tree import DecisionTreeRegressor

X = df[['LSTAT']].values
y = df['MEDV'].values

tree = DecisionTreeRegressor(max_depth=3)
tree.fit(X, y)

sort_idx = X.flatten().argsort()

lin_regplot(X[sort_idx], y[sort_idx], tree)
plt.xlabel('% lower status of the population [LSTAT]')
plt.ylabel('Price in $1000s [MEDV]')
#plt.savefig('images/10_14.png', dpi=300)
plt.show()

结果表明，决策树捕获到了数据中的一般趋势，但是没有捕获到期望预测的连续性和可微性。

8.2随机森林回归

随机森林算法通常比单个决策树具有更好的泛化性能，且有助于降低模型的方差。随机森林算法对于异常值不太敏感，且不需要太多的参数调整。

与使用随机森林分类不同，这里使用MSE准则来生长单个决策树，且预测的目标变量被计算为所有决策树的平均值。

# 使用所用的数据，划分比例为6：4
X = df.iloc[:, :-1].values
y = df['MEDV'].values

X_train, X_test, y_train, y_test = train_test_split(
    X, y, test_size=0.4, random_state=1)

from sklearn.ensemble import RandomForestRegressor

forest = RandomForestRegressor(n_estimators=1000, 
                               criterion='mse', 
                               random_state=1, 
                               n_jobs=-1)
forest.fit(X_train, y_train)
y_train_pred = forest.predict(X_train)
y_test_pred = forest.predict(X_test)

print('MSE train: %.3f, test: %.3f' % (
        mean_squared_error(y_train, y_train_pred),
        mean_squared_error(y_test, y_test_pred)))
print('R^2 train: %.3f, test: %.3f' % (
        r2_score(y_train, y_train_pred),
        r2_score(y_test, y_test_pred)))

MSE train: 1.644, test: 11.085
R^2 train: 0.979, test: 0.877

可以看到，随机森林模型似乎有一定的过拟合，但是在测试集上的$R^2$值为0.87，实际上已经非常高了。

plt.scatter(y_train_pred,  
            y_train_pred - y_train, 
            c='steelblue',
            edgecolor='white',
            marker='o', 
            s=35,
            alpha=0.9,
            label='Training data')
plt.scatter(y_test_pred,  
            y_test_pred - y_test, 
            c='limegreen',
            edgecolor='white',
            marker='s', 
            s=35,
            alpha=0.9,
            label='Test data')

plt.xlabel('Predicted values')
plt.ylabel('Residuals')
plt.legend(loc='upper left')
plt.hlines(y=0, xmin=-10, xmax=50, lw=2, color='black')
plt.xlim([-10, 50])
plt.tight_layout()

#plt.savefig('images/10_15.png', dpi=300)
plt.show()

可以看到，模型对于训练集的拟合依然胜过对测试集的拟合。同时，根据残差图可以看出，残差在零中心点附近的分布似乎不是完全随机的，这表明模型不能

捕捉到所有的特征信息，但与前面各模型的残差图相比，模型性能有了很大的提升。

理想情况下，模型误差应该是随机的或者不可预测的。直观来讲，即预测的误差不应该与特征变量中包含的任何信息有关；

相反，它应该反映出现实世界中分布或模式的随机性。但目前似乎并没有通用的方法来处理残差图中的非随机性，但从建模层面来讲，可以通过特征变换、

超参数调整、选择简单或者更加复杂的模型、去除异常值或者特征提取或构造等方式来改进模型。此外，支持向量回归也常用于非线性问题，文章链接.

9.总结

总结了简单的线性回归分析，对单特征变量和目标变量之间进行了建模。实践了散点图矩阵，查看了数据中的模式和异常；

基于梯度下降的优化方法实现了线性回归模型，然后通过sklearn实现了线性回归建模，以及建立了稳健的回归RANSAC模型；

RANSAC回归可以有效应对异常值；

为了评估模型的性能，计算了模型的MSE和$R^2$指标，此外还对残差图进行了实践操作；

对于线性模型的过拟合，采用了L1和L2正则化，以及两者的折中形式，分别建立了套索回归，岭回归，以及弹性网络模型。

对于非线性问题，应用了多项式回归和随机森林回归，其中一个技巧是对数变换。