线性回归模型是解决回归任务的好起点。
简单线性回归
%matplotlib inline
import matplotlib.pyplot as plt
import seaborn as sns; sns.set()
import numpy as np
rng = np.random.RandomState(1)
x = 10*rng.rand(50)
y = 2*x-5+rng.randn(50)
plt.scatter(x,y);
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);
print("Model slope: ", model.coef_[0])
print("Model intercept:", model.intercept_)
多维度的线性回归模型：
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)
model.intercept_,model.coef_
基函数回归
- 多项式基函数
from sklearn.preprocessing import PolynomialFeatures
x = np.array([2,3,4])
poly = PolynomialFeatures(3,include_bias=False)
poly.fit(x[:,None])
from sklearn.pipeline import make_pipeline
poly_model = make_pipeline(PolynomialFeatures(7),LinearRegression())
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)
- 高斯基函数
from sklearn.base import BaseEstimator, TransformerMixin
class GaussianFeatures(BaseEstimator, TransformerMixin):
"""一维输入均匀分布的高斯特征"""
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):
# 在数据区间中创建N个高斯分布中心
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);
当你在一个Scikit-Learn管道中调用fit方法时，如果管道中包含转换器（如你定义的GaussianFeatures类），那么对于每个转换器，它的fit方法会被自动调用，随后紧接着调用transform方法来转换数据。这个转换后的数据随后传递给管道中的下一个步骤。如果下一步是另一个转换器，这个过程会重复；如果下一步是估计器（如线性回归模型），则只调用估计器的fit方法，因为估计器通常是管道的最后一步。
管道中的这种自动化处理流程是为了简化数据预处理和模型训练的过程。当你调用管道的fit方法时，你实际上是在执行两个主要步骤：
对管道中的每个转换器，依次调用它的fit方法然后是transform方法，以此来训练转换器并转换数据。
对管道中的最后一个步骤（通常是一个估计器），调用它的fit方法，此时使用的数据是经过前面所有转换器处理后的数据。
因此，在你的例子中：gauss_model.fit(x[:, np.newaxis], y)这里gauss_model是一个包含GaussianFeatures转换器的管道。当执行这条fit命令时，实际上会先调用GaussianFeatures的fit方法来找到高斯基函数的中心和宽度，然后调用transform方法将原始数据x转换成基于这些高斯基函数的新特征空间。这个转换后的数据随后被用于管道中的下一个步骤（如线性回归模型）的训练。
这种机制确保了数据可以通过管道中定义的一系列转换步骤被自动处理，最终以适当的形式用于模型训练，而用户无需手动调用每个步骤的fit和transform方法。
案例：自行车流量预测
import pandas as pd
counts = pd.read_csv('SeattleBike-master/FremontHourly.csv', index_col='Date',parse_dates=True)
weather = pd.read_csv('SeattleBike-master/SeaTacWeather.csv', index_col='DATE', parse_dates=True)
daily = counts.resample('d').sum()
daily['Total'] = daily.sum(axis=1)
daily = daily['Total'].to_frame()
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):
"""计算指定日期的白昼时间"""
days = (date - pd.to_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();
# 温度是按照1/10摄氏度统计的，首先转换为摄氏度
weather['TMIN'] /= 10
weather['TMAX'] /= 10
weather['Temp (C)'] = 0.5 * (weather['TMIN'] + weather['TMAX'])
# 降雨量也是按照1/10mm统计的，转化为英寸
weather['PRCP'] /= 254
weather['dry day'] = (weather['PRCP'] == 0).astype(int)
daily = daily.join(weather[['PRCP', 'Temp (C)', 'dry day']])
daily['annual'] = (daily.index - daily.index[0]).days / 365.
print(daily.head())
column_names = ['Mon', 'Tue', 'Wed', 'Thu', 'Fri', 'Sat', 'Sun', 'holiday',
'daylight_hrs', 'PRCP', 'dry day', 'Temp (C)', 'annual']
X = daily[column_names]
y = daily['Total']
model = LinearRegression(fit_intercept=False)
model.fit(X,y)
daily['predicted'] = model.predict(X)
daily[['Total', 'predicted']].plot(alpha=0.5);
params = pd.Series(model.coef_, index=X.columns)
print(params)
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)}))
自举重采样（Bootstrap Resampling）是一种统计方法，用于通过从原始数据中重复抽样（有放回地抽取样本）来估计一个统计量的分布。这种方法特别适用于估计统计量的标准误差或置信区间，特别是在理论分布难以直接计算时。在机器学习和数据科学领域，自举方法可以用来估计模型参数的不确定性。
err = np.std([model.fit(*resample(X, y)).coef_ for i in range(1000)], 0)这行代码实现了自举方法来估计线性回归模型参数的不确定性。具体操作如下：
- 重采样：使用resample(X, y)对原始数据集(X, y)进行重采样，即随机有放回地抽取样本。这意味着某些原始样本可能在重采样的数据集中出现多次，而某些则可能根本不出现。
- 拟合模型：对每个重采样得到的数据集，使用model.fit()方法拟合线性回归模型，并获取模型系数（.coef_）。
- 计算标准差：重复步骤1和2共1000次，对这1000个不同的模型系数集合计算每个系数的标准差。这给出了每个模型参数的不确定性估计，即参数估计的稳定性或变异性。
通过这种方式，err反映了模型参数在不同重采样数据集上的变异程度。如果某个参数的err值较大，意味着该参数的估计值对样本数据较为敏感，其估计的不确定性较高。相反，如果err值较小，则表明参数估计较为稳定，对随机样本抽样的影响较小。
自举重采样是一种强大的非参数统计方法，能够在不依赖于严格假设的前提下，为各种统计量提供有效的不确定性估计，从而在实际应用中非常有用，尤其是在样本大小有限或者理论分布未知的情况下。
在Python中，操作符用于函数调用时，表示将序列（如列表、元组）展开为函数的位置参数。
model.fit(*resample(X, y)).coef_，这里resample函数实际上返回了两个对象（比如两个数组或两个数据帧），操作符的作用就是将这两个返回的对象分别作为model.fit()的第一个参数和第二个参数，而不是将一个包含两个对象的元组作为单一参数传递。
另外，使用随机森林模型可以更准确一些：
from sklearn.ensemble import RandomForestRegressor
model = RandomForestRegressor(n_estimators=100)
# model = LinearRegression(fit_intercept=False)
model.fit(X,y)
daily['predicted'] = model.predict(X)
daily[['Total', 'predicted']].plot(alpha=0.5);
参考：
[1]美 万托布拉斯 (VanderPlas, Jake).Python数据科学手册[M].人民邮电出版社,2018.
