十一、恒定模型、损失和转换
学习成果
- 推导出在 MSE 和 MAE 成本函数下恒定模型的最佳模型参数。
- 评估 MSE 和 MAE 风险之间的差异。
- 理解变量线性化的必要性，并应用图基-莫斯特勒凸图进行转换。
上次，我们介绍了建模过程。我们建立了一个框架，根据一套工作流程，预测目标变量作为我们特征的函数：
- 选择模型 - 我们应该如何表示世界？
- 选择损失函数 - 我们如何量化预测误差？
- 拟合模型 - 我们如何根据我们的数据选择最佳模型参数？
- 评估模型性能 - 我们如何评估这个过程是否产生了一个好模型？
为了说明这个过程，我们推导了简单线性回归（SLR）下均方误差（MSE）作为成本函数的最佳模型参数。SLR 建模过程的摘要如下所示：
在本讲座中，我们将深入探讨步骤 4 - 评估模型性能 - 以 SLR 为例。此外，我们还将通过新模型探索建模过程，继续通过在新模型下找到最佳模型参数来熟悉建模过程，并测试两种不同的损失函数，以了解我们选择的损失如何影响模型设计。稍后，我们将考虑当线性模型不是捕捉数据趋势的最佳选择时会发生什么，以及有哪些解决方案可以创建更好的模型。
11.1 步骤 4：评估 SLR 模型
现在我们已经探讨了（1）选择模型、（2）选择损失函数和（3）拟合模型背后的数学原理，我们还剩下一个最后的问题 - 这个“最佳”拟合模型的预测有多“好”？为了确定这一点，我们可以：
- 可视化数据并计算统计数据：
- 绘制原始数据。
- 计算每一列的均值和标准差。如果我们的预测的均值和标准差接近于原始观察到的y i y_iyi，我们可能会倾向于说我们的模型做得不错。
- （如果我们拟合线性模型）计算相关性r rr。特征和响应变量之间的相关系数的大幅度也可能表明我们的模型做得不错。
- 性能指标：
- 我们可以采用均方根误差（RMSE）。
- 这是均方误差（MSE）的平方根，它是我们一直在最小化以确定最佳模型参数的平均损失。
- RMSE 与y yy的单位相同。
- 较低的 RMSE 表示更“准确”的预测，因为我们在数据中有更低的“平均损失”。
- RMSE = 1 n ∑ i = 1 n ( y i − y ^ i ) 2 \text{RMSE} = \sqrt{\frac{1}{n} \sum_{i=1}^n (y_i - \hat{y}_i)^2}RMSE=n1i=1∑n(yi−y^i)2
- 可视化：
- 查看e i = y i − y i ^ e_i = y_i - \hat{y_i}ei=yi−yi^的残差图，以可视化实际值和预测值之间的差异。良好的残差图不应显示输入/特征x i x_ixi和残差值e i e_iei之间的任何模式。
为了说明这个过程，让我们看看安斯库姆的四重奏。
11.1.1 四个神秘的数据集（安斯库姆的四重奏）
让我们看看四个不同的数据集。
代码
import numpy as np import pandas as pd import matplotlib.pyplot as plt %matplotlib inline import seaborn as sns import itertools from mpl_toolkits.mplot3d import Axes3D
# Big font helper def adjust_fontsize(size=None): SMALL_SIZE = 8 MEDIUM_SIZE = 10 BIGGER_SIZE = 12 if size != None: SMALL_SIZE = MEDIUM_SIZE = BIGGER_SIZE = size plt.rc('font', size=SMALL_SIZE) # controls default text sizes plt.rc('axes', titlesize=SMALL_SIZE) # fontsize of the axes title plt.rc('axes', labelsize=MEDIUM_SIZE) # fontsize of the x and y labels plt.rc('xtick', labelsize=SMALL_SIZE) # fontsize of the tick labels plt.rc('ytick', labelsize=SMALL_SIZE) # fontsize of the tick labels plt.rc('legend', fontsize=SMALL_SIZE) # legend fontsize plt.rc('figure', titlesize=BIGGER_SIZE) # fontsize of the figure title # Helper functions def standard_units(x): return (x - np.mean(x)) / np.std(x) def correlation(x, y): return np.mean(standard_units(x) * standard_units(y)) def slope(x, y): return correlation(x, y) * np.std(y) / np.std(x) def intercept(x, y): return np.mean(y) - slope(x, y)*np.mean(x) def fit_least_squares(x, y): theta_0 = intercept(x, y) theta_1 = slope(x, y) return theta_0, theta_1 def predict(x, theta_0, theta_1): return theta_0 + theta_1*x def compute_mse(y, yhat): return np.mean((y - yhat)**2) plt.style.use('default') # Revert style to default mpl
plt.style.use('default') # Revert style to default mpl NO_VIZ, RESID, RESID_SCATTER = range(3) def least_squares_evaluation(x, y, visualize=NO_VIZ): # statistics print(f"x_mean : {np.mean(x):.2f}, y_mean : {np.mean(y):.2f}") print(f"x_stdev: {np.std(x):.2f}, y_stdev: {np.std(y):.2f}") print(f"r = Correlation(x, y): {correlation(x, y):.3f}") # Performance metrics ahat, bhat = fit_least_squares(x, y) yhat = predict(x, ahat, bhat) print(f"\theta_0: {ahat:.2f}, \theta_1: {bhat:.2f}") print(f"RMSE: {np.sqrt(compute_mse(y, yhat)):.3f}") # visualization fig, ax_resid = None, None if visualize == RESID_SCATTER: fig, axs = plt.subplots(1,2,figsize=(8, 3)) axs[0].scatter(x, y) axs[0].plot(x, yhat) axs[0].set_title("LS fit") ax_resid = axs[1] elif visualize == RESID: fig = plt.figure(figsize=(4, 3)) ax_resid = plt.gca() if ax_resid is not None: ax_resid.scatter(x, y - yhat, color = 'red') ax_resid.plot([4, 14], [0, 0], color = 'black') ax_resid.set_title("Residuals") return fig
# Load in four different datasets: I, II, III, IV x = [10, 8, 13, 9, 11, 14, 6, 4, 12, 7, 5] y1 = [8.04, 6.95, 7.58, 8.81, 8.33, 9.96, 7.24, 4.26, 10.84, 4.82, 5.68] y2 = [9.14, 8.14, 8.74, 8.77, 9.26, 8.10, 6.13, 3.10, 9.13, 7.26, 4.74] y3 = [7.46, 6.77, 12.74, 7.11, 7.81, 8.84, 6.08, 5.39, 8.15, 6.42, 5.73] x4 = [8, 8, 8, 8, 8, 8, 8, 19, 8, 8, 8] y4 = [6.58, 5.76, 7.71, 8.84, 8.47, 7.04, 5.25, 12.50, 5.56, 7.91, 6.89] anscombe = { 'I': pd.DataFrame(list(zip(x, y1)), columns =['x', 'y']), 'II': pd.DataFrame(list(zip(x, y2)), columns =['x', 'y']), 'III': pd.DataFrame(list(zip(x, y3)), columns =['x', 'y']), 'IV': pd.DataFrame(list(zip(x4, y4)), columns =['x', 'y']) } # Plot the scatter plot and line of best fit fig, axs = plt.subplots(2, 2, figsize = (10, 10)) for i, dataset in enumerate(['I', 'II', 'III', 'IV']): ans = anscombe[dataset] x, y = ans['x'], ans['y'] ahat, bhat = fit_least_squares(x, y) yhat = predict(x, ahat, bhat) axs[i//2, i%2].scatter(x, y, alpha=0.6, color='red') # plot the x, y points axs[i//2, i%2].plot(x, yhat) # plot the line of best fit axs[i//2, i%2].set_xlabel(f'$x_{i+1}/details>) axs[i//2, i%2].set_ylabel(f'$y_{i+1}/details>) axs[i//2, i%2].set_title(f"Dataset {dataset}") plt.show()
虽然这四组数据点看起来非常不同，但它们实际上都具有相同的x ˉ \bar xxˉ、y ˉ \bar yyˉ、σ x \sigma_xσx、σ y \sigma_yσy、相关性r rr和 RMSE！如果我们只看这些统计数据，我们可能会倾向于说这些数据集是相似的。
代码
for dataset in ['I', 'II', 'III', 'IV']: print(f">>> Dataset {dataset}:") ans = anscombe[dataset] fig = least_squares_evaluation(ans['x'], ans['y'], visualize = NO_VIZ) print() print()
>>> Dataset I: x_mean : 9.00, y_mean : 7.50 x_stdev: 3.16, y_stdev: 1.94 r = Correlation(x, y): 0.816 heta_0: 3.00, heta_1: 0.50 RMSE: 1.119 >>> Dataset II: x_mean : 9.00, y_mean : 7.50 x_stdev: 3.16, y_stdev: 1.94 r = Correlation(x, y): 0.816 heta_0: 3.00, heta_1: 0.50 RMSE: 1.119 >>> Dataset III: x_mean : 9.00, y_mean : 7.50 x_stdev: 3.16, y_stdev: 1.94 r = Correlation(x, y): 0.816 heta_0: 3.00, heta_1: 0.50 RMSE: 1.118 >>> Dataset IV: x_mean : 9.00, y_mean : 7.50 x_stdev: 3.16, y_stdev: 1.94 r = Correlation(x, y): 0.817 heta_0: 3.00, heta_1: 0.50 RMSE: 1.118
我们可能还希望可视化模型的残差，定义为观察值和预测的y i y_iyi值之间的差异（e i = y i − y ^ i e_i = y_i - \hat{y}_iei=yi−y^i）。这提供了每个预测与真实观察值的“偏差”的高层视图。回想一下，你在Data 8中探讨过这个概念：一个好的回归拟合在其残差图中不应显示出明显的模式。Anscombe 的四重奏的残差图如下所示。请注意，只有第一个图显示出残差大小没有明显模式。这表明 SLR 不是剩下的三组点的最佳模型的指示。
代码
# Residual visualization fig, axs = plt.subplots(2, 2, figsize = (10, 10)) for i, dataset in enumerate(['I', 'II', 'III', 'IV']): ans = anscombe[dataset] x, y = ans['x'], ans['y'] ahat, bhat = fit_least_squares(x, y) yhat = predict(x, ahat, bhat) axs[i//2, i%2].scatter(x, y - yhat, alpha=0.6, color='red') # plot the x, y points axs[i//2, i%2].plot(x, np.zeros_like(x), color='black') # plot the residual line axs[i//2, i%2].set_xlabel(f'$x_{i+1}/details>) axs[i//2, i%2].set_ylabel(f'$e_{i+1}/details>) axs[i//2, i%2].set_title(f"Dataset {dataset} Residuals") plt.show()
11.1.2 预测 vs. 估计
术语预测和估计通常在某种程度上可以互换使用，但它们之间有微妙的区别。估计是使用数据计算模型参数的任务。预测是使用模型预测未见数据的输出的任务。在我们的简单线性回归模型中
y ^ = θ 0 ^ + θ 1 ^ \hat{y} = \hat{\theta_0} + \hat{\theta_1}y^=θ0^+θ1^
我们通过最小化平均损失来估计参数；然后，我们使用这些估计来预测。最小二乘估计是选择最小化 MSE 的参数。
11.2 常数模型 + MSE
现在，我们将从 SLR 模型转换为常数模型，也称为汇总统计。常数模型与我们之前探索过的简单线性回归模型略有不同。常数模型不是从输入的特征变量生成预测，而是始终预测相同的常数数字。这忽略了变量之间的任何关系。例如，假设我们想要预测一家波霸店一天卖出的饮料数量。波霸茶的销售可能取决于一年中的时间、天气、顾客的感觉、学校是否在上课等等，但常数模型忽略了这些因素，而更倾向于一个更简单的模型。换句话说，常数模型采用了一个简化的假设。
它也是一个参数化的统计模型：
y ^ i = θ 0 \hat{y}_i = \theta_0y^i=θ0
θ 0 \theta_0θ0是常数模型的参数，就像θ 0 \theta_0θ0和θ 1 \theta_1θ1是 SLR 中的参数一样。由于我们的参数θ 0 \theta_0θ0是一维的（θ 0 ∈ R \theta_0 \in \mathbb{R}θ0∈R），我们现在的模型没有输入，将始终预测y ^ i = θ 0 \hat{y}_i = \theta_0y^i=θ0。
11.2.1 推导最优的θ 0 \theta_0θ0
我们现在的任务是确定什么值的θ 0 \theta_0θ0最能代表最佳模型 - 换句话说，每次猜测什么数字可以在我们的数据上获得最低可能的平均损失？
像以前一样，我们将使用均方误差（MSE）。回想一下，MSE 是数据D = { y 1 , y 2 , . . . , y n } D = \{y_1, y_2, ..., y_n\}D={y1,y2,...,yn}上的平均平方损失（L2 损失）。
R ( θ ) = 1 n ∑ i = 1 n ( y i − y i ^ ) 2 R(\theta) = \frac{1}{n}\sum^{n}_{i=1} (y_i - \hat{y_i})^2R(θ)=n1i=1∑n(yi−yi^)2
我们的建模过程现在看起来像这样：
- 选择模型：常数模型
- 选择损失函数：L2 损失
- 拟合模型
- 评估模型性能
UCB Data100：数据科学的原理和技巧：第十一章到第十二章（2）https://developer.aliyun.com/article/1427175
