5.6 線性回歸
譯者：飛龍
協(xié)議：CC BY-NC-SA 4.0
譯文沒有得到原作者授權(quán)，不保證與原文的意思嚴(yán)格一致。
就像樸素貝葉斯（之前在樸素貝葉斯分類中討論）是分類任務(wù)的一個(gè)很好的起點(diǎn)，線性回歸模型是回歸任務(wù)的一個(gè)很好的起點(diǎn)。 這些模型受歡迎，因?yàn)樗鼈兛梢钥焖贁M合，并且非?？山忉尅?你可能熟悉線性回歸模型的最簡單形式（即使用直線擬合數(shù)據(jù)），但是可以擴(kuò)展這些模型，來建模更復(fù)雜的數(shù)據(jù)行為。
在本節(jié)中，在這個(gè)眾所周知問題背后，我們將從數(shù)學(xué)的快速直觀的了解開始，然后再看看如何將線性模型推廣到數(shù)據(jù)中更復(fù)雜的模式。
我們以標(biāo)準(zhǔn)導(dǎo)入來開始：
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns; sns.set()
import numpy as np
簡單線性回歸
我們以最熟悉的線性回歸開始，它是一個(gè)擬合數(shù)據(jù)的直線。擬合直線是這樣：
y = ax + b
其中a是眾所周知的斜率，b是眾所周知的截距。
考慮下面的數(shù)據(jù)，它點(diǎn)綴在直線周圍，斜率為 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估計(jì)其來擬合這個(gè)直線，并且構(gòu)造出最佳擬合直線。
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);
數(shù)據(jù)的斜率和截距包含在模型的擬合參數(shù)中，其中 Scikit-Learn 總是以尾部的下劃線標(biāo)記。 這里的相關(guān)參數(shù)是coef_和intercept_：
print("Model slope: ", model.coef_[0])
print("Model intercept:", model.intercept_)
Model slope: 2.02720881036
Model intercept: -4.99857708555
我們看到結(jié)果非常接近輸入，這是我們希望的。
然而，線性回歸估計(jì)器比這更加強(qiáng)大，除了簡單的直線擬合之外，它還可以處理這種形式的多維線性模型。
y = a0 + a1x1 + a2x2 + ...
其中有多個(gè)x值。 在幾何學(xué)上，這類似于使用平面擬合三維點(diǎn)，或使用超平面擬合更高維度的點(diǎn)。
這種回歸的多維本質(zhì)使得它們更難以可視化，但是我們可以通過使用 NumPy 的矩陣乘法運(yùn)算符，構(gòu)建一些示例數(shù)據(jù)來查看其中的一個(gè)擬合：
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. ]
這里數(shù)據(jù)y由三個(gè)隨機(jī)x值構(gòu)成，線性回歸恢復(fù)了用于構(gòu)建數(shù)據(jù)的系數(shù)。
以這種方式，我們可以使用單個(gè)LinearRegression估計(jì)器來將數(shù)據(jù)擬合為直線，平面或超平面。 這種方法似乎仍然限制于變量之間的嚴(yán)格線性關(guān)系，但事實(shí)證明，我們也可以使其寬松。
基函數(shù)回歸
用于將線性回歸適配變量之間的非線性關(guān)系的一個(gè)技巧是，根據(jù)基函數(shù)來轉(zhuǎn)換數(shù)據(jù)。在超參數(shù)和模型驗(yàn)證和特征工程中使用的PolynomialRegression流水線中，我們已經(jīng)看到了其中的一個(gè)版本。這個(gè)想法是使用我們的多維線性模型：
y = a0 + a1x1 + a2x2 + ...
并從我們的一維輸入x中構(gòu)造x1、x2、x3，以及其他。也就是，我們讓xn = fn(x)，其中fn是轉(zhuǎn)換數(shù)據(jù)的函數(shù)。
例如，如果fn(x) = x ** n，我們的模型就變成了多項(xiàng)式回歸：
y = a0 + a1x + a2x^2 + a3x^3 + ...
要注意這仍然是個(gè)線性模型 -- 線性是指，參數(shù)an永遠(yuǎn)不會互相相乘或者相除。我們所做的是選取我們的一維x值并投影到更高維度上，以便線性擬合可以擬合x和y的更復(fù)雜關(guān)系。
多項(xiàng)式基函數(shù)
多項(xiàng)式投影足夠?qū)嵱?，它?nèi)建于 Scikit-Learn，使用PolynomialFeatures轉(zhuǎn)換器。
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.]])
我們在這里看到，通過取每個(gè)值的指數(shù)，轉(zhuǎn)換器將我們的一維數(shù)組轉(zhuǎn)換為三維數(shù)組。 然后這種新的更高維度的數(shù)據(jù)表示，可以用于線性回歸。
我們在特征工程中看到，實(shí)現(xiàn)這一目標(biāo)的最簡單的方法是使用流水線。我們以這種方式制作一個(gè) 7 階多項(xiàng)式模型：
from sklearn.pipeline import make_pipeline
poly_model = make_pipeline(PolynomialFeatures(7),
LinearRegression())
通過這種轉(zhuǎn)換，我們可以使用線性模型來擬合x和y之間更復(fù)雜的關(guān)系。 例如，這里是一個(gè)帶有噪音的正弦波：
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 階的多項(xiàng)式基函數(shù)，我們的線性模型可以提供這個(gè)非線性數(shù)據(jù)的優(yōu)秀的擬合。
高斯基函數(shù)
當(dāng)然，也存在其他的基函數(shù)。 例如，一個(gè)有用的模式是擬合一個(gè)模型，它不是多項(xiàng)式基數(shù)的和，而是高斯基數(shù)的和。 結(jié)果可能如下圖所示：
繪圖中的陰影區(qū)域是縮放的基函數(shù)，并且當(dāng)它們相加時(shí)，它們通過數(shù)據(jù)再現(xiàn)平滑曲線。 這些高斯基函數(shù)不內(nèi)置在 Scikit-Learn 中，但是我們可以編寫一個(gè)自定義的轉(zhuǎn)換器來創(chuàng)建它們，如下圖所示（Scikit-Learn 轉(zhuǎn)換器實(shí)現(xiàn)為 Python 類；閱讀 Scikit-Learn 的源代碼，是了解如何創(chuàng)建它們的好方法）：
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);
我們把這個(gè)例子放在這里，只是為了說明多項(xiàng)式基函數(shù)沒有任何魔法：如果你對數(shù)據(jù)的生成過程有一些直覺，它使你認(rèn)為一個(gè)或另一個(gè)基可能是適當(dāng)?shù)?，你就可以使用它們?/p>
正則化
將基函數(shù)引入到我們的線性回歸中，使得模型更加靈活，但也可以很快導(dǎo)致過擬合（參考在超參數(shù)和模型驗(yàn)證和特征工程中的討論）。 例如，如果我們選擇太多的高斯基函數(shù)，我們最終得到的結(jié)果看起來不太好：
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);
將數(shù)據(jù)投影到 30 維的基上，該模型具有太多的靈活性，并且在由數(shù)據(jù)約束的位置之間達(dá)到極值。 如果我們繪制高斯基相對于它們的位置的系數(shù)，我們可以看到這個(gè)原因：
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)
該圖的底部圖像顯示了基函數(shù)在每個(gè)位置的幅度。 當(dāng)基函數(shù)重疊時(shí)，這是典型的過擬合行為：相鄰基函數(shù)的系數(shù)相互排斥并相互抵消。 我們知道這樣的行為是有問題的，如果我們可以通過懲罰模型參數(shù)的較大值，來限制模型中的這種尖峰。 這種懲罰被稱為正則化，有幾種形式。
嶺回歸（L2 正則化）
也許最常見的正則化形式被稱為嶺回歸或 L2 正則化，有時(shí)也稱為 Tikhonov 正則化。 這通過懲罰模型系數(shù)的平方和（第二范數(shù)）來實(shí)現(xiàn)；在這種情況下，模型擬合的懲罰將是：
其中α是控制懲罰強(qiáng)度的自由參數(shù)。 這種類型的懲罰模型是構(gòu)建在 Scikit-Learn 中的Ridge估計(jì)器：
from sklearn.linear_model import Ridge
model = make_pipeline(GaussianFeatures(30), Ridge(alpha=0.1))
basis_plot(model, title='Ridge Regression')
α參數(shù)基本上是控制所得模型復(fù)雜度的旋鈕。 在極限α→0中，結(jié)果恢復(fù)標(biāo)準(zhǔn)線性回歸; 在極限α→∞中，所有模型響應(yīng)都將被抑制。 特別是嶺回歸的一個(gè)優(yōu)點(diǎn)是，它的計(jì)算很高效 -- 比原始線性回歸模型，幾乎不需要更多的計(jì)算成本。
套索回歸（L1 正則化）
另一種非常常見的正則化類型稱為套索，懲罰回歸系數(shù)的絕對值（第一范數(shù)）之和：
雖然這在概念上非常類似于嶺回歸，但是結(jié)果可能會令人驚訝：例如，由于幾何原因，套索回歸可能的情況下傾向于稀疏模型：即，它優(yōu)先將模型系數(shù)設(shè)置為恰好為零。
我們可以復(fù)制嶺回歸的圖像，但使用 L1 歸一化系數(shù)，看到這種行為：
from sklearn.linear_model import Lasso
model = make_pipeline(GaussianFeatures(30), Lasso(alpha=0.001))
basis_plot(model, title='Lasso Regression')
利用套索回歸的乘法，大多數(shù)系數(shù)恰好為零，功能行為由可用基函數(shù)的一小部分建模。 與嶺正則化一樣，α參數(shù)調(diào)整懲罰的強(qiáng)度，并且應(yīng)通過例如交叉驗(yàn)證來確定（參考超參數(shù)和模型驗(yàn)證和特征工程中的討論）。
示例：預(yù)測自行車流量
例如，我們來看看我們是否可以根據(jù)天氣，季節(jié)和其他因素，來預(yù)測西雅圖 Fremont 大橋的自行車流量。 我們已經(jīng)在使用時(shí)間序列中，看到這些數(shù)據(jù)。
在本節(jié)中，我們將把自行車數(shù)據(jù)與另一個(gè)數(shù)據(jù)連接到一起，并嘗試確定天氣和季節(jié)因素（溫度，降水和日光時(shí)間）在多大程度上影響這條路上的自行車流量。 幸運(yùn)的是，NOAA 提供他們的氣象站日常數(shù)據(jù)（我使用站點(diǎn)號碼 USW00024233），我們可以輕松地使用 Pandas 連接兩個(gè)數(shù)據(jù)源。 我們將執(zhí)行一個(gè)簡單的線性回歸，將天氣和其他信息與自行車計(jì)數(shù)相關(guān)聯(lián)，以便估計(jì)這些參數(shù)中的任何一個(gè)的變化，如何影響特定日期的人數(shù)。
特別地，這是一個(gè)例子，說明如何將 Scikit-Learn 的工具用于統(tǒng)計(jì)建模框架，其中假定模型的參數(shù)具有可解釋的含義。 如前所述，這不是機(jī)器學(xué)習(xí)中的標(biāo)準(zhǔn)方法，但是對于某些模型，可以這么解釋。
我們開始加載兩個(gè)數(shù)據(jù)集，按日期索引：
# !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)
下面我們會計(jì)算總共的自行車日常流量，并將其放到自己的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)
與之類似，我們可能預(yù)期，騎車人在假期表現(xiàn)不同。讓我們?yōu)榇艘蔡砑右粋€(gè)指標(biāo)。
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)
我們也猜測，日光的小時(shí)數(shù)也應(yīng)該騎車人數(shù)。讓我們使用標(biāo)準(zhǔn)的天文計(jì)算來添加這個(gè)信息。
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)
我們還可以向數(shù)據(jù)中加入平均溫度和總降水。 除了降水的英寸之外，我們還添加一個(gè)標(biāo)志，指示一天是否干燥（降雨量為零）：
# 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']])
最后，讓我們添加一個(gè)從第 1 天起增加的計(jì)數(shù)器，并且測量已經(jīng)過去了幾年。 這將讓我們度量，在每日的流量中，任何觀測到的年度增長或下降：
daily['annual'] = (daily.index - daily.index[0]).days / 365.
現(xiàn)在我們的數(shù)據(jù)有序了，我們可以看一看：
daily.head()
最后，我們可以在視覺上，比較總共的和預(yù)測的自行車流量。
daily[['Total', 'predicted']].plot(alpha=0.5);
很明顯，我們錯(cuò)過了一些關(guān)鍵特征，特別是在夏季。 我們的特征還不完整（即，人們不僅僅根據(jù)這些，決定是否騎車去上班），或者有一些非線性關(guān)系，我們沒有考慮到（例如，也許人們在高和低溫度下騎行較少）。 然而，我們的粗略近似足以提供一些見解，我們可以看一下線性模型的系數(shù)，來估計(jì)每個(gè)特征對每日自行車數(shù)量的貢獻(xiàn)：
params = pd.Series(model.coef_, index=X.columns)
params
Mon 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
沒有一些不確定性的度量，這些數(shù)字難以解釋。 我們可以使用數(shù)據(jù)的引導(dǎo)重采樣，來快速計(jì)算這些不確定性：
from sklearn.utils import resample
np.random.seed(1)
err = np.std([model.fit(*resample(X, y)).coef_
for i in range(1000)], 0)
估計(jì)了這些誤差，讓我們再次看看結(jié)果：
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
我們首先看到，每周的基線的趨勢相對穩(wěn)定：平日比周末和假日有更多的騎車人。我們看到，每增加一小時(shí)的日光，就多出129 ± 9個(gè)騎車人; 每升高 1 攝氏度的溫度，就多出65 ± 4人；干燥的日子意味著平均有548 ± 33個(gè)騎車人，每一英寸的降水意味著665 ± ??62個(gè)人將自行車放在家里。一旦考慮到所有這些影響，我們每年都會適度增加27 ± 18新的日常騎車人。
我們的模型幾乎肯定缺少一些相關(guān)信息。例如，這個(gè)模型不能解釋非線性效應(yīng)（如降水和低溫的影響）以及每個(gè)變量內(nèi)的非線性趨勢（如在非常冷和非常熱的溫度下的騎行傾向）。此外，我們已經(jīng)拋棄了一些更細(xì)致的信息（如雨天的早上和下午之間的差異），我們忽略了天數(shù)之間的相關(guān)性（例如星期二下雨可能影響周三的數(shù)值，或連續(xù)下雨后的意想不到的陽光燦爛的日子的效果）。這些都是潛在的有趣效果，如果你愿意，你現(xiàn)在有了開始探索它們的工具！
