Python 数据科学手册 5.6 线性回归
5.6 线性回归
原文：In Depth: Linear Regression
译者：飞龙
译文没有得到原作者授权，不保证与原文的意思严格一致。
就像朴素贝叶斯（之前在朴素贝叶斯分类中讨论）是分类任务的一个很好的起点，线性回归模型是回归任务的一个很好的起点。 这些模型受欢迎，因为它们可以快速拟合，并且非常可解释。 你可能熟悉线性回归模型的最简单形式（即使用直线拟合数据），但是可以扩展这些模型，来建模更复杂的数据行为。
在本节中，在这个众所周知问题背后，我们将从数学的快速直观的了解开始，然后再看看如何将线性模型推广到数据中更复杂的模式。
我们以标准导入来开始：
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns; sns.set()
import numpy as np简单线性回归
我们以最熟悉的线性回归开始，它是一个拟合数据的直线。拟合直线是这样：
y = ax + b其中a是众所周知的斜率，b是众所周知的截距。
考虑下面的数据，它点缀在直线周围，斜率为 2，截距为 -5。
rng = np.random.RandomState(1)
x = 10 * rng.rand(50)
y = 2 * x - 5 + rng.randn(50)
plt.scatter(x, y);
我们可以使用 Scikit-Learn 的LinearRegression估计其来拟合这个直线，并且构造出最佳拟合直线。
from sklearn.linear_model import LinearRegression
model = LinearRegression(fit_intercept=True)
model.fit(x[:, np.newaxis], y)
xfit = np.linspace(0, 10, 1000)
yfit = model.predict(xfit[:, np.newaxis])
plt.scatter(x, y)
plt.plot(xfit, yfit);
数据的斜率和截距包含在模型的拟合参数中，其中 Scikit-Learn 总是以尾部的下划线标记。 这里的相关参数是coef_和intercept_：
print("Model slope: ", model.coef_[0])
print("Model intercept:", model.intercept_)Model slope: 2.02720881036
Model intercept: -4.99857708555我们看到结果非常接近输入，这是我们希望的。
然而，线性回归估计器比这更加强大，除了简单的直线拟合之外，它还可以处理这种形式的多维线性模型。
y = a0 + a1x1 + a2x2 + ...其中有多个x值。 在几何学上，这类似于使用平面拟合三维点，或使用超平面拟合更高维度的点。
这种回归的多维本质使得它们更难以可视化，但是我们可以通过使用 NumPy 的矩阵乘法运算符，构建一些示例数据来查看其中的一个拟合：
rng = np.random.RandomState(1)
X = 10 * rng.rand(100, 3)
y = 0.5 + np.dot(X, [1.5, -2., 1.])
model.fit(X, y)
print(model.intercept_)
print(model.coef_)0.5
[ 1.5 -2. 1. ]这里数据y由三个随机x值构成，线性回归恢复了用于构建数据的系数。
以这种方式，我们可以使用单个LinearRegression估计器来将数据拟合为直线，平面或超平面。 这种方法似乎仍然限制于变量之间的严格线性关系，但事实证明，我们也可以使其宽松。
基函数回归
用于将线性回归适配变量之间的非线性关系的一个技巧是，根据基函数来转换数据。在超参数和模型验证和特征工程中使用的PolynomialRegression流水线中，我们已经看到了其中的一个版本。这个想法是使用我们的多维线性模型：
y = a0 + a1x1 + a2x2 + ...并从我们的一维输入x中构造x1、x2、x3，以及其他。也就是，我们让xn = fn(x)，其中fn是转换数据的函数。
例如，如果fn(x) = x ** n，我们的模型就变成了多项式回归：
y = a0 + a1x + a2x^2 + a3x^3 + ...要注意这仍然是个线性模型 – 线性是指，参数an永远不会互相相乘或者相除。我们所做的是选取我们的一维x值并投影到更高维度上，以便线性拟合可以拟合x和y的更复杂关系。
多项式基函数
多项式投影足够实用，它内建于 Scikit-Learn，使用PolynomialFeatures转换器。
from sklearn.preprocessing import PolynomialFeatures
x = np.array([2, 3, 4])
poly = PolynomialFeatures(3, include_bias=False)
poly.fit_transform(x[:, None])array([[ 2., 4., 8.],
[ 3., 9., 27.],
[ 4., 16., 64.]])我们在这里看到，通过取每个值的指数，转换器将我们的一维数组转换为三维数组。 然后这种新的更高维度的数据表示，可以用于线性回归。
我们在特征工程中看到，实现这一目标的最简单的方法是使用流水线。我们以这种方式制作一个 7 阶多项式模型：
from sklearn.pipeline import make_pipeline
poly_model = make_pipeline(PolynomialFeatures(7),
LinearRegression())通过这种转换，我们可以使用线性模型来拟合x和y之间更复杂的关系。 例如，这里是一个带有噪音的正弦波：
rng = np.random.RandomState(1)
x = 10 * rng.rand(50)
y = np.sin(x) + 0.1 * rng.randn(50)
poly_model.fit(x[:, np.newaxis], y)
yfit = poly_model.predict(xfit[:, np.newaxis])
plt.scatter(x, y)
plt.plot(xfit, yfit);
使用 7 阶的多项式基函数，我们的线性模型可以提供这个非线性数据的优秀的拟合。
高斯基函数
当然，也存在其他的基函数。 例如，一个有用的模式是拟合一个模型，它不是多项式基数的和，而是高斯基数的和。 结果可能如下图所示：
绘图中的阴影区域是缩放的基函数，并且当它们相加时，它们通过数据再现平滑曲线。 这些高斯基函数不内置在 Scikit-Learn 中，但是我们可以编写一个自定义的转换器来创建它们，如下图所示（Scikit-Learn 转换器实现为 Python 类；阅读 Scikit-Learn 的源代码，是了解如何创建它们的好方法）：
from sklearn.base import BaseEstimator, TransformerMixin
class GaussianFeatures(BaseEstimator, TransformerMixin):
"""Uniformly spaced Gaussian features for one-dimensional input"""
def __init__(self, N, width_factor=2.0):
self.N = N
self.width_factor = width_factor
@staticmethod
def _gauss_basis(x, y, width, axis=None):
arg = (x - y) / width
return np.exp(-0.5 * np.sum(arg ** 2, axis))
def fit(self, X, y=None):
# create N centers spread along the data range
self.centers_ = np.linspace(X.min(), X.max(), self.N)
self.width_ = self.width_factor * (self.centers_[1] - self.centers_[0])
return self
def transform(self, X):
return self._gauss_basis(X[:, :, np.newaxis], self.centers_,
self.width_, axis=1)
gauss_model = make_pipeline(GaussianFeatures(20),
LinearRegression())
gauss_model.fit(x[:, np.newaxis], y)
yfit = gauss_model.predict(xfit[:, np.newaxis])
plt.scatter(x, y)
plt.plot(xfit, yfit)
plt.xlim(0, 10);
我们把这个例子放在这里，只是为了说明多项式基函数没有任何魔法：如果你对数据的生成过程有一些直觉，它使你认为一个或另一个基可能是适当的，你就可以使用它们。
正则化
将基函数引入到我们的线性回归中，使得模型更加灵活，但也可以很快导致过拟合（参考在超参数和模型验证和特征工程中的讨论）。 例如，如果我们选择太多的高斯基函数，我们最终得到的结果看起来不太好：
model = make_pipeline(GaussianFeatures(30),
LinearRegression())
model.fit(x[:, np.newaxis], y)
plt.scatter(x, y)
plt.plot(xfit, model.predict(xfit[:, np.newaxis]))
plt.xlim(0, 10)
plt.ylim(-1.5, 1.5);
将数据投影到 30 维的基上，该模型具有太多的灵活性，并且在由数据约束的位置之间达到极值。 如果我们绘制高斯基相对于它们的位置的系数，我们可以看到这个原因：
def basis_plot(model, title=None):
fig, ax = plt.subplots(2, sharex=True)
model.fit(x[:, np.newaxis], y)
ax[0].scatter(x, y)
ax[0].plot(xfit, model.predict(xfit[:, np.newaxis]))
ax[0].set(xlabel='x', ylabel='y', ylim=(-1.5, 1.5))
if title:
ax[0].set_title(title)
ax[1].plot(model.steps[0][1].centers_,
model.steps[1][1].coef_)
ax[1].set(xlabel='basis location',
ylabel='coefficient',
xlim=(0, 10))
model = make_pipeline(GaussianFeatures(30), LinearRegression())
basis_plot(model)
该图的底部图像显示了基函数在每个位置的幅度。 当基函数重叠时，这是典型的过拟合行为：相邻基函数的系数相互排斥并相互抵消。 我们知道这样的行为是有问题的，如果我们可以通过惩罚模型参数的较大值，来限制模型中的这种尖峰。 这种惩罚被称为正则化，有几种形式。
岭回归（L2 正则化）
也许最常见的正则化形式被称为岭回归或 L2 正则化，有时也称为 Tikhonov 正则化。 这通过惩罚模型系数的平方和（第二范数）来实现；在这种情况下，模型拟合的惩罚将是：
其中α是控制惩罚强度的自由参数。 这种类型的惩罚模型是构建在 Scikit-Learn 中的Ridge估计器：
from sklearn.linear_model import Ridge
model = make_pipeline(GaussianFeatures(30), Ridge(alpha=0.1))
basis_plot(model, title='Ridge Regression')
α参数基本上是控制所得模型复杂度的旋钮。 在极限α→0中，结果恢复标准线性回归; 在极限α→∞中，所有模型响应都将被抑制。 特别是岭回归的一个优点是，它的计算很高效 – 比原始线性回归模型，几乎不需要更多的计算成本。
套索回归（L1 正则化）
另一种非常常见的正则化类型称为套索，惩罚回归系数的绝对值（第一范数）之和：
虽然这在概念上非常类似于岭回归，但是结果可能会令人惊讶：例如，由于几何原因，套索回归可能的情况下倾向于稀疏模型：即，它优先将模型系数设置为恰好为零。
我们可以复制岭回归的图像，但使用 L1 归一化系数，看到这种行为：
from sklearn.linear_model import Lasso
model = make_pipeline(GaussianFeatures(30), Lasso(alpha=0.001))
basis_plot(model, title='Lasso Regression')
利用套索回归的乘法，大多数系数恰好为零，功能行为由可用基函数的一小部分建模。 与岭正则化一样，α参数调整惩罚的强度，并且应通过例如交叉验证来确定（参考超参数和模型验证和特征工程中的讨论）。
示例：预测自行车流量
例如，我们来看看我们是否可以根据天气，季节和其他因素，来预测西雅图 Fremont 大桥的自行车流量。 我们已经在使用时间序列中，看到这些数据。
在本节中，我们将把自行车数据与另一个数据连接到一起，并尝试确定天气和季节因素（温度，降水和日光时间）在多大程度上影响这条路上的自行车流量。 幸运的是，NOAA 提供他们的气象站日常数据（我使用站点号码 USW00024233），我们可以轻松地使用 Pandas 连接两个数据源。 我们将执行一个简单的线性回归，将天气和其他信息与自行车计数相关联，以便估计这些参数中的任何一个的变化，如何影响特定日期的人数。
特别地，这是一个例子，说明如何将 Scikit-Learn 的工具用于统计建模框架，其中假定模型的参数具有可解释的含义。 如前所述，这不是机器学习中的标准方法，但是对于某些模型，可以这么解释。
我们开始加载两个数据集，按日期索引：
# !curl -o FremontBridge.csv https://data.seattle.gov/api/views/65db-xm6k/rows.csv?accessType=DOWNLOAD
import pandas as pd
counts = pd.read_csv('FremontBridge.csv', index_col='Date', parse_dates=True)
weather = pd.read_csv('data/BicycleWeather.csv', index_col='DATE', parse_dates=True)下面我们会计算总共的自行车日常流量，并将其放到自己的DataFrame中：
daily = counts.resample('d').sum()
daily['Total'] = daily.sum(axis=1)
daily = daily[['Total']] # remove other columns我们以前看到，使用模式通常每天都有所不同; 我们通过添加表示星期几的二元列，来解释它：
days = ['Mon', 'Tue', 'Wed', 'Thu', 'Fri', 'Sat', 'Sun']
for i in range(7):
daily[days[i]] = (daily.index.dayofweek == i).astype(float)与之类似，我们可能预期，骑车人在假期表现不同。让我们为此也添加一个指标。
from pandas.tseries.holiday import USFederalHolidayCalendar
cal = USFederalHolidayCalendar()
holidays = cal.holidays('2012', '2016')
daily = daily.join(pd.Series(1, index=holidays, name='holiday'))
daily['holiday'].fillna(0, inplace=True)我们也猜测，日光的小时数也应该骑车人数。让我们使用标准的天文计算来添加这个信息。
def hours_of_daylight(date, axis=23.44, latitude=47.61):
"""Compute the hours of daylight for the given date"""
days = (date - pd.datetime(2000, 12, 21)).days
m = (1. - np.tan(np.radians(latitude))
* np.tan(np.radians(axis) * np.cos(days * 2 * np.pi / 365.25)))
return 24. * np.degrees(np.arccos(1 - np.clip(m, 0, 2))) / 180.
daily['daylight_hrs'] = list(map(hours_of_daylight, daily.index))
daily[['daylight_hrs']].plot()
plt.ylim(8, 17)
# (8, 17)
我们还可以向数据中加入平均温度和总降水。 除了降水的英寸之外，我们还添加一个标志，指示一天是否干燥（降雨量为零）：
# temperatures are in 1/10 deg C; convert to C
weather['TMIN'] /= 10
weather['TMAX'] /= 10
weather['Temp (C)'] = 0.5 * (weather['TMIN'] + weather['TMAX'])
# precip is in 1/10 mm; convert to inches
weather['PRCP'] /= 254
weather['dry day'] = (weather['PRCP'] == 0).astype(int)
daily = daily.join(weather[['PRCP', 'Temp (C)', 'dry day']])最后，让我们添加一个从第 1 天起增加的计数器，并且测量已经过去了几年。 这将让我们度量，在每日的流量中，任何观测到的年度增长或下降：
daily['annual'] = (daily.index - daily.index[0]).days / 365.现在我们的数据有序了，我们可以看一看：
daily.head()
最后，我们可以在视觉上，比较总共的和预测的自行车流量。
daily[['Total', 'predicted']].plot(alpha=0.5);
很明显，我们错过了一些关键特征，特别是在夏季。 我们的特征还不完整（即，人们不仅仅根据这些，决定是否骑车去上班），或者有一些非线性关系，我们没有考虑到（例如，也许人们在高和低温度下骑行较少）。 然而，我们的粗略近似足以提供一些见解，我们可以看一下线性模型的系数，来估计每个特征对每日自行车数量的贡献：
params = pd.Series(model.coef_, index=X.columns)
paramsMon 504.882756
Tue 610.233936
Wed 592.673642
Thu 482.358115
Fri 177.980345
Sat -1103.301710
Sun -1133.567246
holiday -1187.401381
daylight_hrs 128.851511
PRCP -664.834882
dry day 547.698592
Temp (C) 65.162791
annual 26.942713
dtype: float64没有一些不确定性的度量，这些数字难以解释。 我们可以使用数据的引导重采样，来快速计算这些不确定性：
from sklearn.utils import resample
np.random.seed(1)
err = np.std([model.fit(*resample(X, y)).coef_
for i in range(1000)], 0)估计了这些误差，让我们再次看看结果：
print(pd.DataFrame({'effect': params.round(0),
'error': err.round(0)})) effect error
Mon 505.0 86.0
Tue 610.0 83.0
Wed 593.0 83.0
Thu 482.0 85.0
Fri 178.0 81.0
Sat -1103.0 80.0
Sun -1134.0 83.0
holiday -1187.0 163.0
daylight_hrs 129.0 9.0
PRCP -665.0 62.0
dry day 548.0 33.0
Temp (C) 65.0 4.0
annual 27.0 18.0我们首先看到，每周的基线的趋势相对稳定：平日比周末和假日有更多的骑车人。我们看到，每增加一小时的日光，就多出129 ± 9个骑车人; 每升高 1 摄氏度的温度，就多出65 ± 4人；干燥的日子意味着平均有548 ± 33个骑车人，每一英寸的降水意味着665 ± 62个人将自行车放在家里。一旦考虑到所有这些影响，我们每年都会适度增加27 ± 18新的日常骑车人。
我们的模型几乎肯定缺少一些相关信息。例如，这个模型不能解释非线性效应（如降水和低温的影响）以及每个变量内的非线性趋势（如在非常冷和非常热的温度下的骑行倾向）。此外，我们已经抛弃了一些更细致的信息（如雨天的早上和下午之间的差异），我们忽略了天数之间的相关性（例如星期二下雨可能影响周三的数值，或连续下雨后的意想不到的阳光灿烂的日子的效果）。这些都是潜在的有趣效果，如果你愿意，你现在有了开始探索它们的工具！
