5.6 線性回歸
譯者：飛龍
協(xié)議：CC BY-NC-SA 4.0
譯文沒有得到原作者授權(quán)，不保證與原文的意思嚴格一致。
就像樸素貝葉斯（之前在樸素貝葉斯分類中討論）是分類任務(wù)的一個很好的起點璃饱，線性回歸模型是回歸任務(wù)的一個很好的起點。 這些模型受歡迎，因為它們可以快速擬合滓技，并且非常可解釋棚潦。 你可能熟悉線性回歸模型的最簡單形式（即使用直線擬合數(shù)據(jù)）令漂，但是可以擴展這些模型，來建模更復(fù)雜的數(shù)據(jù)行為丸边。
在本節(jié)中叠必，在這個眾所周知問題背后，我們將從數(shù)學(xué)的快速直觀的了解開始原环，然后再看看如何將線性模型推廣到數(shù)據(jù)中更復(fù)雜的模式挠唆。
我們以標準導(dǎo)入來開始：
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns; sns.set()
import numpy as np
簡單線性回歸
我們以最熟悉的線性回歸開始，它是一個擬合數(shù)據(jù)的直線嘱吗。擬合直線是這樣：
y = ax + b
其中a是眾所周知的斜率玄组，b是眾所周知的截距。
考慮下面的數(shù)據(jù)谒麦，它點綴在直線周圍俄讹，斜率為 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估計其來擬合這個直線患膛，并且構(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 總是以尾部的下劃線標記踪蹬。 這里的相關(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é)果非常接近輸入，這是我們希望的臣咖。
然而跃捣，線性回歸估計器比這更加強大，除了簡單的直線擬合之外夺蛇，它還可以處理這種形式的多維線性模型疚漆。
y = a0 + a1x1 + a2x2 + ...
其中有多個x值。 在幾何學(xué)上，這類似于使用平面擬合三維點娶聘，或使用超平面擬合更高維度的點闻镶。
這種回歸的多維本質(zhì)使得它們更難以可視化，但是我們可以通過使用 NumPy 的矩陣乘法運算符丸升，構(gòu)建一些示例數(shù)據(jù)來查看其中的一個擬合：
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由三個隨機x值構(gòu)成铆农，線性回歸恢復(fù)了用于構(gòu)建數(shù)據(jù)的系數(shù)。
以這種方式发钝，我們可以使用單個LinearRegression估計器來將數(shù)據(jù)擬合為直線顿涣，平面或超平面。 這種方法似乎仍然限制于變量之間的嚴格線性關(guān)系酝豪，但事實證明，我們也可以使其寬松精堕。
基函數(shù)回歸
用于將線性回歸適配變量之間的非線性關(guān)系的一個技巧是孵淘，根據(jù)基函數(shù)來轉(zhuǎn)換數(shù)據(jù)。在超參數(shù)和模型驗證和特征工程中使用的PolynomialRegression流水線中歹篓，我們已經(jīng)看到了其中的一個版本瘫证。這個想法是使用我們的多維線性模型：
y = a0 + a1x1 + a2x2 + ...
并從我們的一維輸入x中構(gòu)造x1、x2庄撮、x3背捌，以及其他。也就是洞斯，我們讓xn = fn(x)毡庆，其中fn是轉(zhuǎn)換數(shù)據(jù)的函數(shù)。
例如烙如，如果fn(x) = x ** n么抗，我們的模型就變成了多項式回歸：
y = a0 + a1x + a2x^2 + a3x^3 + ...
要注意這仍然是個線性模型 -- 線性是指，參數(shù)an永遠不會互相相乘或者相除亚铁。我們所做的是選取我們的一維x值并投影到更高維度上蝇刀，以便線性擬合可以擬合x和y的更復(fù)雜關(guān)系。
多項式基函數(shù)
多項式投影足夠?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.]])
我們在這里看到然爆，通過取每個值的指數(shù)站粟，轉(zhuǎn)換器將我們的一維數(shù)組轉(zhuǎn)換為三維數(shù)組。 然后這種新的更高維度的數(shù)據(jù)表示施蜜，可以用于線性回歸卒蘸。
我們在特征工程中看到，實現(xiàn)這一目標的最簡單的方法是使用流水線。我們以這種方式制作一個 7 階多項式模型：
from sklearn.pipeline import make_pipeline
poly_model = make_pipeline(PolynomialFeatures(7),
LinearRegression())
通過這種轉(zhuǎn)換缸沃，我們可以使用線性模型來擬合x和y之間更復(fù)雜的關(guān)系恰起。 例如，這里是一個帶有噪音的正弦波：
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 階的多項式基函數(shù)趾牧，我們的線性模型可以提供這個非線性數(shù)據(jù)的優(yōu)秀的擬合检盼。
高斯基函數(shù)
當(dāng)然，也存在其他的基函數(shù)翘单。 例如吨枉，一個有用的模式是擬合一個模型，它不是多項式基數(shù)的和哄芜，而是高斯基數(shù)的和貌亭。 結(jié)果可能如下圖所示：
繪圖中的陰影區(qū)域是縮放的基函數(shù)，并且當(dāng)它們相加時认臊，它們通過數(shù)據(jù)再現(xiàn)平滑曲線圃庭。 這些高斯基函數(shù)不內(nèi)置在 Scikit-Learn 中，但是我們可以編寫一個自定義的轉(zhuǎn)換器來創(chuàng)建它們失晴，如下圖所示（Scikit-Learn 轉(zhuǎn)換器實現(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);
我們把這個例子放在這里涂屁，只是為了說明多項式基函數(shù)沒有任何魔法：如果你對數(shù)據(jù)的生成過程有一些直覺书在，它使你認為一個或另一個基可能是適當(dāng)?shù)模憔涂梢允褂盟鼈儭?/p>
正則化
將基函數(shù)引入到我們的線性回歸中拆又，使得模型更加靈活儒旬，但也可以很快導(dǎo)致過擬合（參考在超參數(shù)和模型驗證和特征工程中的討論）。 例如遏乔，如果我們選擇太多的高斯基函數(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ù)約束的位置之間達到極值凉翻。 如果我們繪制高斯基相對于它們的位置的系數(shù)，我們可以看到這個原因：
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ù)在每個位置的幅度捻激。 當(dāng)基函數(shù)重疊時制轰，這是典型的過擬合行為：相鄰基函數(shù)的系數(shù)相互排斥并相互抵消。 我們知道這樣的行為是有問題的胞谭，如果我們可以通過懲罰模型參數(shù)的較大值垃杖，來限制模型中的這種尖峰。 這種懲罰被稱為正則化丈屹，有幾種形式调俘。
嶺回歸（L2 正則化）
也許最常見的正則化形式被稱為嶺回歸或 L2 正則化伶棒，有時也稱為 Tikhonov 正則化。 這通過懲罰模型系數(shù)的平方和（第二范數(shù)）來實現(xiàn)彩库；在這種情況下肤无，模型擬合的懲罰將是：
其中α是控制懲罰強度的自由參數(shù)。 這種類型的懲罰模型是構(gòu)建在 Scikit-Learn 中的Ridge估計器：
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ù)標準線性回歸; 在極限α→∞中，所有模型響應(yīng)都將被抑制眯搭。 特別是嶺回歸的一個優(yōu)點是窥翩，它的計算很高效 -- 比原始線性回歸模型，幾乎不需要更多的計算成本鳞仙。
套索回歸（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)整懲罰的強度，并且應(yīng)通過例如交叉驗證來確定（參考超參數(shù)和模型驗證和特征工程中的討論）掂墓。
示例：預(yù)測自行車流量
例如谦纱，我們來看看我們是否可以根據(jù)天氣，季節(jié)和其他因素君编，來預(yù)測西雅圖 Fremont 大橋的自行車流量跨嘉。 我們已經(jīng)在使用時間序列中，看到這些數(shù)據(jù)吃嘿。
在本節(jié)中祠乃，我們將把自行車數(shù)據(jù)與另一個數(shù)據(jù)連接到一起，并嘗試確定天氣和季節(jié)因素（溫度兑燥，降水和日光時間）在多大程度上影響這條路上的自行車流量亮瓷。 幸運的是，NOAA 提供他們的氣象站日常數(shù)據(jù)（我使用站點號碼 USW00024233）降瞳，我們可以輕松地使用 Pandas 連接兩個數(shù)據(jù)源嘱支。 我們將執(zhí)行一個簡單的線性回歸，將天氣和其他信息與自行車計數(shù)相關(guān)聯(lián)，以便估計這些參數(shù)中的任何一個的變化除师，如何影響特定日期的人數(shù)沛膳。
特別地，這是一個例子馍盟，說明如何將 Scikit-Learn 的工具用于統(tǒng)計建挠谥茫框架，其中假定模型的參數(shù)具有可解釋的含義贞岭。 如前所述八毯，這不是機器學(xué)習(xí)中的標準方法，但是對于某些模型瞄桨，可以這么解釋话速。
我們開始加載兩個數(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)
下面我們會計算總共的自行車日常流量芯侥，并將其放到自己的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)榇艘蔡砑右粋€指標唉工。
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ù)也應(yīng)該騎車人數(shù)。讓我們使用標準的天文計算來添加這個信息淋硝。
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ù)中加入平均溫度和總降水雹熬。 除了降水的英寸之外，我們還添加一個標志谣膳，指示一天是否干燥（降雨量為零）：
# 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 天起增加的計數(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);
很明顯捂襟，我們錯過了一些關(guān)鍵特征，特別是在夏季欢峰。 我們的特征還不完整（即葬荷，人們不僅僅根據(jù)這些涨共，決定是否騎車去上班），或者有一些非線性關(guān)系宠漩，我們沒有考慮到（例如举反，也許人們在高和低溫度下騎行較少）。 然而扒吁，我們的粗略近似足以提供一些見解火鼻，我們可以看一下線性模型的系數(shù)，來估計每個特征對每日自行車數(shù)量的貢獻：
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)重采樣，來快速計算這些不確定性：
from sklearn.utils import resample
np.random.seed(1)
err = np.std([model.fit(*resample(X, y)).coef_
for i in range(1000)], 0)
估計了這些誤差盼铁，讓我們再次看看結(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)定：平日比周末和假日有更多的騎車人。我們看到饶火，每增加一小時的日光鹏控，就多出129 ± 9個騎車人; 每升高 1 攝氏度的溫度，就多出65 ± 4人肤寝；干燥的日子意味著平均有548 ± 33個騎車人当辐，每一英寸的降水意味著665 ± ??62個人將自行車放在家里。一旦考慮到所有這些影響鲤看，我們每年都會適度增加27 ± 18新的日常騎車人瀑构。
我們的模型幾乎肯定缺少一些相關(guān)信息。例如刨摩，這個模型不能解釋非線性效應(yīng)（如降水和低溫的影響）以及每個變量內(nèi)的非線性趨勢（如在非常冷和非常熱的溫度下的騎行傾向）。此外世吨，我們已經(jīng)拋棄了一些更細致的信息（如雨天的早上和下午之間的差異）澡刹，我們忽略了天數(shù)之間的相關(guān)性（例如星期二下雨可能影響周三的數(shù)值，或連續(xù)下雨后的意想不到的陽光燦爛的日子的效果）耘婚。這些都是潛在的有趣效果罢浇，如果你愿意，你現(xiàn)在有了開始探索它們的工具沐祷！
