Python之ML–回归分析预测连续型目标变量
监督学习的另一个分支：回归分析(regression analysis).回归模型(regression model)可用于连续型目标变量的预测分析
主要知识点如下：- 数据集的探索与可视化
- 实现线性回归模型的不同方法
- 训练可处理异常值的回归模型
- 回归模型的评估及常见问题
- 基于非线性数据拟合回归模型
一.简单线性回归模型
简单(单变量)线性回归的目标是：通过模型来描述某一特征(解释变量x)与连续输出(目标变量y)之间的关系.当只有一个解释变量时,线性模型的函数定义如下：
其中,权值w0为函数在y轴上的截距,w1为解释变量的系数.我们的目标是通过学习得到线性方程的这两个权值,并用它们描述解释变量与目标变量之间的关系 我们可以将线性回归模型扩展为多个解释变量.即为所谓的多元线性回归(multiple linear regression)
from IPython.display import Image
其中,w0为x0=1时在y轴上的截距
二.波士顿房屋数据集
波士顿房屋数据集(Housing Dataset),可以在UCI机器学习数据库中下载：
数据集中506个样本的特征描述如下：
- CRIM:房屋所在镇的犯罪率
- ZN:用地面积远大于25000平方英尺的住宅所占比例
- INDUS:房屋所在镇无零售业务区域所占比例
- CHAS:与查尔斯河有关的虚拟变量(如果房屋位于河边则为1,否则为0)
- NOX:一氧化碳浓度
- RM:每处寓所的平均房间数
- AGE:业主自住房屋中,建于1940年之前的房屋所占比例
- DIS:房屋距离波士顿五大就业中心的加权距离
- RAD:距离房屋最近公路入口编号
- TAX:每一万美元全额财产税金额
- PTRATIO:房屋所在镇的师生比
- B:计算公式为1000(Bk-0.63)^2,其中Bk为房屋所在镇非裔美籍人口所占比例
- LSTAT:弱势群体人口所占比例
- MEDV:业主自住房屋的平均价格
先从UCI库中把Housing Dataset导入到pandas的DataFrame中：
import pandas as pddf = pd.read_csv('https://archive.ics.uci.edu/ml/machine-learning-databases/' 'housing/housing.data', header=None, sep='\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 |
df = pd.read_csv('https://raw.githubusercontent.com/rasbt/python-machine-learning-book/master/code/datasets/housing/housing.data', header=None, sep='\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 |
首先,借助散点图矩阵,我们以可视化的方法汇总显示各不同特征两两之间的关系.为了绘制散点图矩阵,我们需要用到seaborn库中的pairplot函数,它是在matplotlib基础上绘制统计图像的python库
import matplotlib.pyplot as pltimport seaborn as snssns.set(style='whitegrid', context='notebook')cols = ['LSTAT', 'INDUS', 'NOX', 'RM', 'MEDV']sns.pairplot(df[cols], size=2.5)plt.tight_layout()plt.show()
通过此散点图矩阵,我们可以快速了解数据是如何分布的,以及其中是否包含异常值.我们可直观看出RM和房屋价格MEDV(第5列和第4行)之间存在线性关系.从MEDV直方图中可以发现：MEDV看似呈正态分布,但包含几个异常值
使用Numpy的corrcoef函数计算前面散点图矩阵中5个特征间的相关系数矩阵,并使用seaborn的heatmap函数绘制其对应的热度图
import numpy as npcm = np.corrcoef(df[cols].values.T)sns.set(font_scale=1.5)hm = sns.heatmap(cm, cbar=True, annot=True, square=True, fmt='.2f', annot_kws={ 'size': 15}, yticklabels=cols, xticklabels=cols)plt.show() 
为了拟合线性回归模型,我们主要关注那些更目标变量MEDV高度相关的特征.观察前面的相关系数矩阵,可以发现MEDV与变量LSTAT的相关性最大(-0.74)。前面的散点图矩阵显示LSTAT和MEDV之间存在明显的非线性关系.另一方面,RM和MEDV间的相关性也较高(0.70)
导入seaborn库后,会覆盖当前python会话中matplotlib默认的图像显示方式.如果不希望使用seaborn风格的设置,可以通过如下命令重设为matplotlib的风格：
sns.reset_orig()
二.基于最小二乘法构建线性回归模型
通过最小二乘法(Ordinary Least Squares,OLS)估计回归曲线的参数,使得回归曲线到样本点垂直距离(残差或误差)的平方和最小
1.通过梯度下降计算回归参数
Adaline中的代价函数就是误差平方和(Sum of Squared Error,SSE),它等同于我们定义的OLS代价函数.本质上,OLS线性回归可以理解为无单位阶跃函数的Adaline,这样我们得到的是连续型的输出值,而不是以-1和1来代表的类标.为了说明两者的相似程度,使用前面实现的GD方法,并且移除其单位阶跃函数来实现我们第一个线性回归模型：
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']].valuesy = df['MEDV'].values
from sklearn.preprocessing import StandardScalersc_x = StandardScaler()sc_y = StandardScaler()X_std = sc_x.fit_transform(X)y_std = sc_y.fit_transform(y[:, np.newaxis]).flatten()
lr = LinearRegressionGD()lr.fit(X_std, y_std)
<__main__.LinearRegressionGD at 0x278bdee8908>
我们将再次绘制迭代次数对应的代价函数的值,以检查线性回归是否收敛
plt.plot(range(1, lr.n_iter+1), lr.cost_)plt.ylabel('SSE')plt.xlabel('Epoch')plt.tight_layout()plt.show() 
我们使用lin_regplot函数绘制房间数与房屋价格之间的关系图
def lin_regplot(X,y,model): plt.scatter(X,y,c='blue') plt.plot(X,model.predict(X),color='red') return None
lin_regplot(X_std,y_std,lr)plt.xlabel('Average number of rows [RM] (Standardized)')plt.ylabel('Price in $1000\'s [MEDV] (standardized)')plt.show() 
为了将预测价格缩放到以1000美元为价格单位的坐标轴中,我们使用了StandardScaler的inverse_transform方法
num_rooms_std=sc_x.transform(np.array([[5.0]]))price_std=lr.predict(num_rooms_std)print("Price in $1000's:%.3f"%sc_y.inverse_transform(price_std)) Price in $1000's:10.840
print('Slope:%.3f'%lr.w_[1])print('Intercept:.%.3f'%lr.w_[0]) Slope:0.695Intercept:.-0.000
2.使用scikit-learn估计回归模型的系数
from sklearn.linear_model import LinearRegressionslr=LinearRegression()slr.fit(X,y)print('Slope:%.3f'%slr.coef_[0])print('Intercept:%.3f'%slr.intercept_) Slope:9.102Intercept:-34.671
lin_regplot(X,y,slr)plt.xlabel('Average number of rows [RM] (Standardized)')plt.ylabel('Price in $1000\'s [MEDV] (standardized)')plt.show() 
上述代码绘制出了训练数据和拟合模型的图像,总体效果与GD算法实现的模型是一致的
四.使用RANSAC拟合高鲁棒性回归模型
作为清除异常值的一种高鲁棒性回归方法,在此我们将学习随机抽样一致性(RANdom SAmple Consensus,RANSAC)算法,使用数据的一个子集(即所谓的内点,inlier)来进行回归模型的拟合
我们将RANSAC算法的工作流程总结如下：
- 从数据集中随机抽取样本构建内点集合来拟合模型
- 使用剩余数据对上一步得到的模型进行测试,并将落在预定公差范围内的样本点增至内点集合中
- 使用全部的内点集合数据再次进行模型的拟合
- 使用内点集合来估计模型的误差
- 如果模型性能达到了用户设定的特定阈值或者迭代达到了预定次数,则算法终止,否则跳转到第1步
from sklearn.linear_model import RANSACRegressorransac=RANSACRegressor(LinearRegression(), max_trials=100, min_samples=50, loss='absolute_loss', residual_threshold=5.0, random_state=0)ransac.fit(X,y)
RANSACRegressor(base_estimator=LinearRegression(copy_X=True, fit_intercept=True, n_jobs=None, normalize=False), is_data_valid=None, is_model_valid=None, loss='absolute_loss', max_skips=inf, max_trials=100, min_samples=50, random_state=0, residual_threshold=5.0, stop_n_inliers=inf, stop_probability=0.99, stop_score=inf)
我们将RANSACRegression的最大迭代次数设定为100,设定参数min_samples=50,即随机抽取作为内点的最小样本数量设定为50.将residual_threshold参数的值设定为5.0,这样只有与拟合曲线垂直距离小于5个单位长度的样本点才能加入到内点集合,此设置在特定数据集上表现良好
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='blue', marker='o', label='Inliers')plt.scatter(X[outlier_mask], y[outlier_mask], c='lightgreen', marker='s', label='Outliers')plt.plot(line_X, line_y_ransac, color='red') plt.xlabel('Average number of rooms [RM]')plt.ylabel('Price in $1000\'s [MEDV]')plt.legend(loc='upper left')plt.tight_layout()plt.show() 
显示模型的斜率和截距
print('Slope:%.3f'%ransac.estimator_.coef_[0])print('Intercept:.%.3f'%ransac.estimator_.intercept_) Slope:10.735Intercept:.-44.089
五.线性回归模型性能的评估
使用数据集中的所有变量训练多元回归模型
from sklearn.model_selection import train_test_splitX=df.iloc[:,:-1].valuesy=df['MEDV'].valuesX_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)
我们绘制得到残差图,其中通过将预测结果减去对应目标变量真实值,便可获得残差的值
plt.scatter(y_train_pred,y_train_pred-y_train,c='blue',marker='o',label='Training data')plt.scatter(y_test_pred,y_test_pred-y_test,c='lightgreen',marker='s',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='red')plt.xlim([-10,50])plt.show() 
from sklearn.metrics import mean_squared_errorprint('MSE train:%.3f,test:%.3f'%(mean_squared_error(y_train,y_train_pred),mean_squared_error(y_test,y_test_pred))) MSE train:19.958,test:27.196
from sklearn.metrics import r2_scoreprint('R^2 train:%.3f,test:%.3f'%(r2_score(y_train,y_train_pred),r2_score(y_test,y_test_pred))) R^2 train:0.765,test:0.673
六.线性回归模型的曲线化–多项式回归
使用scikit-learn中的PolynomialFeatures转换类在只含一个解释变量的简单回归问题中加入二次项(d=2),并且将多项式回归与线性回归进行线性拟合
1.增加一个二次多项式项
from sklearn.preprocessing import PolynomialFeaturesX=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])lr=LinearRegression()pr=LinearRegression()quadratic=PolynomialFeatures(degree=2)X_quad=quadratic.fit_transform(X)
2.拟合一个用于对比的简单线性回归模型
lr.fit(X,y)X_fit=np.arange(250,600,10)[:,np.newaxis]y_lin_fit=lr.predict(X_fit)
3.使用经过转换后的特征针对多项式回归拟合一个多元线性回归模型
pr.fit(X_quad, y)y_quad_fit = pr.predict(quadratic.fit_transform(X_fit))# plot resultsplt.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.legend(loc='upper left')plt.tight_layout()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.330Training R^2 linear: 0.832, quadratic: 0.982
1.房屋数据集中的非线性关系建模
我们将使用二次和三次多项式对房屋价格和LSTAT(弱势群体人口所占比例)之间的关系进行建模,并与线性拟合进行对比
X = df[['LSTAT']].valuesy = df['MEDV'].valuesregr = LinearRegression()# create quadratic featuresquadratic = PolynomialFeatures(degree=2)cubic = PolynomialFeatures(degree=3)X_quad = quadratic.fit_transform(X)X_cubic = cubic.fit_transform(X)# fit featuresX_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))# plot resultsplt.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 $1000\'s [MEDV]')plt.legend(loc='upper right')plt.tight_layout()plt.show() 
X = df[['LSTAT']].valuesy = df['MEDV'].values# transform featuresX_log = np.log(X)y_sqrt = np.sqrt(y)# fit featuresX_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))# plot resultsplt.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 \; \$1000\'s [MEDV]}$')plt.legend(loc='lower left')plt.tight_layout()plt.show() 
2.使用随机森林处处理非线性关系
随机森林(random forest)是多棵决策树(decision tree)的集合,与先前介绍的线性回归和多项式回归不同,它可以被理解为分段线性函数的集成.换句话说,通过决策树算法,我们把输入空间细分更小的区域以更好地管理
决策树回归
决策树算法的一个优点在于：如果我们处理的是非线性数据,无需对其进行特征转换
使用DecisionTreeRegressor类对MEDV和LSTAT两个变量之间的非线性关系进行建模
from sklearn.tree import DecisionTreeRegressorX = df[['LSTAT']].valuesy = df['MEDV'].valuestree = 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 $1000\'s [MEDV]')plt.show() 
随机森林回归
由于随机性有助于降低模型的方差,与单棵决策树相比,随机森林通常具有更好的泛化性能.随机森林的另一个优势在于：它对数据集中的异常值不敏感,且无需过多的参数调优.随机森林中唯一需要我们通过实验来获得的参数就是待集成决策树的数量.唯一区别在于：随机森林回归使用MSE作为单棵决策树生成的标准,同时所有决策树预测值的平均数作为预测目标变量的值
我们使用房屋数据集中的所有特征来拟合一个随机森林回归模型,其中60%的样本用于模型的训练,剩余的40%用来对模型进行评估
X = df.iloc[:, :-1].valuesy = df['MEDV'].valuesX_train, X_test, y_train, y_test = train_test_split( X, y, test_size=0.4, random_state=1)
from sklearn.ensemble import RandomForestRegressorforest = 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.641, test: 11.056R^2 train: 0.979, test: 0.878
plt.scatter(y_train_pred, y_train_pred - y_train, c='black', marker='o', s=35, alpha=0.5, label='Training data')plt.scatter(y_test_pred, y_test_pred - y_test, c='lightgreen', marker='s', s=35, alpha=0.7, 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='red')plt.xlim([-10, 50])plt.tight_layout()plt.show() 