十一、恒定模型、损失和转换

原文：Constant Model, Loss, and Transformations

译者：飞龙

协议：CC BY-NC-SA 4.0

学习成果

• 推导出在 MSE 和 MAE 成本函数下恒定模型的最佳模型参数。

• 评估 MSE 和 MAE 风险之间的差异。

• 理解变量线性化的必要性，并应用图基-莫斯特勒凸图进行转换。

上次，我们介绍了建模过程。我们建立了一个框架，根据一套工作流程，预测目标变量作为我们特征的函数：

1. 选择模型 - 我们应该如何表示世界？

2. 选择损失函数 - 我们如何量化预测误差？

3. 拟合模型 - 我们如何根据我们的数据选择最佳模型参数？

4. 评估模型性能 - 我们如何评估这个过程是否产生了一个好模型？

为了说明这个过程，我们推导了简单线性回归（SLR）下均方误差（MSE）作为成本函数的最佳模型参数。SLR 建模过程的摘要如下所示：

在本讲座中，我们将深入探讨步骤 4 - 评估模型性能 - 以 SLR 为例。此外，我们还将通过新模型探索建模过程，继续通过在新模型下找到最佳模型参数来熟悉建模过程，并测试两种不同的损失函数，以了解我们选择的损失如何影响模型设计。稍后，我们将考虑当线性模型不是捕捉数据趋势的最佳选择时会发生什么，以及有哪些解决方案可以创建更好的模型。

11.1 步骤 4：评估 SLR 模型

现在我们已经探讨了（1）选择模型、（2）选择损失函数和（3）拟合模型背后的数学原理，我们还剩下一个最后的问题 - 这个“最佳”拟合模型的预测有多“好”？为了确定这一点，我们可以：

1. 可视化数据并计算统计数据：

   • 绘制原始数据。

   • 计算每一列的均值和标准差。如果我们的预测的均值和标准差接近于原始观察到的$y_i$，我们可能会倾向于说我们的模型做得不错。

   • （如果我们拟合线性模型）计算相关性$r$。特征和响应变量之间的相关系数的大幅度也可能表明我们的模型做得不错。

2. 性能指标：

   • 我们可以采用均方根误差（RMSE）。

     • 这是均方误差（MSE）的平方根，它是我们一直在最小化以确定最佳模型参数的平均损失。

     • RMSE 与$y$的单位相同。

     • 较低的 RMSE 表示更“准确”的预测，因为我们在数据中有更低的“平均损失”。

   $\text{RMSE} = \sqrt{\frac{1}{n} \sum_{i=1}^n (y_i - \hat{y}_i)^2}$

3. 可视化：

   • 查看$e_i = y_i - \hat{y_i}$的残差图，以可视化实际值和预测值之间的差异。良好的残差图不应显示输入/特征$x_i$和残差值$e_i$之间的任何模式。

为了说明这个过程，让我们看看安斯库姆的四重奏。

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()

虽然这四组数据点看起来非常不同，但它们实际上都具有相同的$\bar x$、$\bar y$、$\sigma_x$、$\sigma_y$、相关性$r$和 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$值之间的差异（$e_i = y_i - \hat{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. 估计

术语预测和估计通常在某种程度上可以互换使用，但它们之间有微妙的区别。估计是使用数据计算模型参数的任务。预测是使用模型预测未见数据的输出的任务。在我们的简单线性回归模型中

$\hat{y} = \hat{\theta_0} + \hat{\theta_1}$

我们通过最小化平均损失来估计参数；然后，我们使用这些估计来预测。最小二乘估计是选择最小化 MSE 的参数。

11.2 常数模型 + MSE

现在，我们将从 SLR 模型转换为常数模型，也称为汇总统计。常数模型与我们之前探索过的简单线性回归模型略有不同。常数模型不是从输入的特征变量生成预测，而是始终预测相同的常数数字。这忽略了变量之间的任何关系。例如，假设我们想要预测一家波霸店一天卖出的饮料数量。波霸茶的销售可能取决于一年中的时间、天气、顾客的感觉、学校是否在上课等等，但常数模型忽略了这些因素，而更倾向于一个更简单的模型。换句话说，常数模型采用了一个简化的假设。

它也是一个参数化的统计模型：

$\hat{y}_i = \theta_0$

$\theta_0$是常数模型的参数，就像$\theta_0$和$\theta_1$是 SLR 中的参数一样。由于我们的参数$\theta_0$是一维的（$\theta_0 \in \mathbb{R}$），我们现在的模型没有输入，将始终预测$\hat{y}_i = \theta_0$。

11.2.1 推导最优的$\theta_0$

我们现在的任务是确定什么值的$\theta_0$最能代表最佳模型 - 换句话说，每次猜测什么数字可以在我们的数据上获得最低可能的平均损失？

像以前一样，我们将使用均方误差（MSE）。回想一下，MSE 是数据$D = \{y_1, y_2, ..., y_n\}$上的平均平方损失（L2 损失）。

$R(\theta) = \frac{1}{n}\sum^{n}_{i=1} (y_i - \hat{y_i})^2$

我们的建模过程现在看起来像这样：

1. 选择模型：常数模型

2. 选择损失函数：L2 损失

3. 拟合模型

4. 评估模型性能

给定常数模型$\hat{y}_i = \theta_0$，我们可以将 MSE 方程重写为

$R(\theta) = \frac{1}{n}\sum^{n}_{i=1} (y_i - \theta_0)^2$

我们可以通过找到最优的$\theta_0$来拟合模型，从而最小化 MSE，使用微积分方法。

1. 对$\theta_0$求导

\begin{align} \frac{d}{d\theta_0}\text{R}(\theta) & = \frac{d}{d\theta_0}\frac{1}{n}\sum^{n}_{i=1} (y_i - \theta_0)^2 \\ &= {n}\sum^{n}_{i=1} \frac{d}{d\theta_0} (y_i - \theta_0)^2 \quad \quad \text{求和的导数是导数的和} \\ &= {n}\sum^{n}_{i=1} 2 (y_i - \theta_0) (-1) \quad \quad \text{链式法则} \\ &= {\frac{-2}{n}}\sum^{n}_{i=1} (y_i - \theta_0) \quad \quad \text{简单的常数} \end{align}

1. 等于 0 $0 = {\frac{-2}{n}}\sum^{n}_{i=1} (y_i - \theta_0)$

2. 解出 $\theta_0$

\begin{align} 0 &= {\frac{-2}{n}}\sum^{n}_{i=1} (y_i - \theta_0) \\ &= \sum^{n}_{i=1} (y_i - \theta_0) \quad \quad \text{两边同时除以} \frac{-2}{n} \\ &= \sum^{n}_{i=1} y_i - \sum^{n}_{i=1} \theta_0 \quad \quad \text{分开求和} \\ &= \sum^{n}_{i=1} y_i - n * \theta_0 \quad \quad \text{c + c + … + c = nc} \\ n * \theta_0 &= \sum^{n}_{i=1} y_i \\ \theta_0 &= \frac{1}{n} \sum^{n}_{i=1} y_i \\ \theta_0 &= \bar{y} \end{align}

让我们花点时间解释一下这个结果。 $\hat{\theta} = \bar{y}$ 是常数模型 + MSE 的最佳参数。无论你有什么样的数据样本，它都是成立的，并且它提供了一些正式的推理，解释了为什么均值是如此常见的摘要统计量。

我们的最佳模型参数是使成本函数最小化的参数值。成本函数的最小值可以表示为：

$R(\hat{\theta}) = \min_{\theta} R(\theta)$

用简单的英语重新陈述上面的内容：当成本函数以最佳参数作为输入时，我们正在查看成本函数的值。这个最佳模型参数 $\hat{\theta}$ 是使成本 $R$ 最小化的 $\theta$ 的值。

对于建模目的，我们更关心成本的最小值 $R(\hat{\theta})$，而不是导致这种最低平均损失的 * $\theta$ 的值。换句话说，我们关心找到最佳参数值，使得：

$\hat{\theta} = \underset{\theta}{\operatorname{\arg\min}}\:R(\theta)$

也就是说，我们想要找到使成本函数最小化的参数 $\theta$。

11.2.2 比较两个不同的模型，都使用 MSE 进行拟合

现在我们已经探讨了带有 L2 损失的常数模型，我们可以将其与上一讲学到的 SLR 模型进行比较。考虑下面的数据集，其中包含嘴海牛的年龄和长度信息。假设我们想要预测嘴海牛的年龄：

	常数模型	简单线性回归
模型	$\hat{y} = \theta_0$	$\hat{y} = \theta_0 + \theta1 x$
数据	年龄样本 $D = \{y_1, y_2, ..., y_m\}$	年龄样本 $D = \{(x_1, y_1), (x_2, y_2), ..., (x_n, y_n)\}$
维度	$\hat{\theta_0}$ 是 1-D	$\hat{\theta} = [\hat{\theta_0}, \hat{\theta_1}]$ 是 2-D
损失曲面	2-D	3-D
损失模型	$\hat{R}(\theta) = \frac{1}{n}\sum^{n}_{i=1} (y_i - \theta_0)^2$	$\hat{R}(\theta) = \frac{1}{n}\sum^{n}_{i=1} (y_i - (\theta_0 + \theta_1 x))^2$
RMSE	7.72	4.31
可视化预测	地毯图	散点图

（注意我们的 SLR 散点图的点在视觉上并不是一个很好的线性拟合。我们会回到这个问题）。

生成图形和模型的代码如下，但我们不会深入讨论。

代码

dugongs = pd.read_csv("data/dugongs.csv")
data_constant = dugongs["Age"]
data_linear = dugongs[["Length", "Age"]]

# Constant Model + MSE
plt.style.use('default') # Revert style to default mpl
adjust_fontsize(size=16)
%matplotlib inline

def mse_constant(theta, data):
 return np.mean(np.array([(y_obs - theta) ** 2 for y_obs in data]), axis=0)

thetas = np.linspace(-20, 42, 1000)
l2_loss_thetas = mse_constant(thetas, data_constant)

# Plotting the loss surface
plt.plot(thetas, l2_loss_thetas)
plt.xlabel(r'$\theta_0/details>)
plt.ylabel(r'MSE')

# Optimal point
thetahat = np.mean(data_constant)
plt.scatter([thetahat], [mse_constant(thetahat, data_constant)], s=50, label = r"$\hat{\theta}_0$")
plt.legend();
# plt.show()

代码

# SLR + MSE
def mse_linear(theta_0, theta_1, data_linear):
 data_x, data_y = data_linear.iloc[:,0], data_linear.iloc[:,1]
 return np.mean(np.array([(y - (theta_0+theta_1*x)) ** 2 for x, y in zip(data_x, data_y)]), axis=0)

# plotting the loss surface
theta_0_values = np.linspace(-80, 20, 80)
theta_1_values = np.linspace(-10, 30, 80)
mse_values = np.array([[mse_linear(x,y,data_linear) for x in theta_0_values] for y in theta_1_values])

# Optimal point
data_x, data_y = data_linear.iloc[:, 0], data_linear.iloc[:, 1]
theta_1_hat = np.corrcoef(data_x, data_y)[0, 1] * np.std(data_y) / np.std(data_x)
theta_0_hat = np.mean(data_y) - theta_1_hat * np.mean(data_x)

# Create the 3D plot
fig = plt.figure(figsize=(7, 5))
ax = fig.add_subplot(111, projection='3d')

X, Y = np.meshgrid(theta_0_values, theta_1_values)
surf = ax.plot_surface(X, Y, mse_values, cmap='viridis', alpha=0.6)  # Use alpha to make it slightly transparent

# Scatter point using matplotlib
sc = ax.scatter([theta_0_hat], [theta_1_hat], [mse_linear(theta_0_hat, theta_1_hat, data_linear)],
 marker='o', color='red', s=100, label='theta hat')

# Create a colorbar
cbar = fig.colorbar(surf, ax=ax, shrink=0.5, aspect=10)
cbar.set_label('Cost Value')

ax.set_title('MSE for different $\\theta_0, \\theta_1/details>)
ax.set_xlabel('$\\theta_0/details>)
ax.set_ylabel('$\\theta_1/details>) 
ax.set_zlabel('MSE')

# plt.show()

Text(0.5, 0, 'MSE')

代码

# Predictions
yobs = data_linear["Age"]      # The true observations y
xs = data_linear["Length"]     # Needed for linear predictions
n = len(yobs)                  # Predictions

yhats_constant = [thetahat for i in range(n)]    # Not used, but food for thought
yhats_linear = [theta_0_hat + theta_1_hat * x for x in xs]

# Constant Model Rug Plot
# In case we're in a weird style state
sns.set_theme()
adjust_fontsize(size=16)
%matplotlib inline

fig = plt.figure(figsize=(8, 1.5))
sns.rugplot(yobs, height=0.25, lw=2) ;
plt.axvline(thetahat, color='red', lw=4, label=r"$\hat{\theta}_0$");
plt.legend()
plt.yticks([]);
# plt.show()

代码

# SLR model scatter plot 
# In case we're in a weird style state
sns.set_theme()
adjust_fontsize(size=16)
%matplotlib inline

sns.scatterplot(x=xs, y=yobs)
plt.plot(xs, yhats_linear, color='red', lw=4);
# plt.savefig('dugong_line.png', bbox_inches = 'tight');
# plt.show()

解释 RMSE（均方根误差）：* 常数误差高于线性误差。

因此，* 常数模型比线性模型更差（至少对于这个度量）。

11.3 常数模型 + MAE

我们现在看到，改变用于预测的模型会导致最佳模型参数的结果大不相同。如果我们改变模型评估中使用的损失函数会发生什么？

这一次，我们将考虑具有 L1（绝对损失）作为损失函数的常数模型。这意味着平均损失将被表示为平均绝对误差（MAE）。

1. 选择模型：常数模型

2. 选择损失函数：L1 损失

3. 拟合模型

4. 评估模型性能

11.3.1 求解最优 $\theta_0$

回想一下，MAE 是数据 $D = \{y_1, y_2, ..., y_m\}$ 上的平均绝对损失（L1 损失）。

$\hat{R}(\theta) = \frac{1}{n}\sum^{n}_{i=1} |y_i - \hat{y_i}|$

给定常数模型 $\hat{y} = \theta_0$，我们可以将 MAE 写成：

$\hat{R}(\theta) = \frac{1}{n}\sum^{n}_{i=1} |y_i - \theta_0|$

为了拟合模型，我们通过微积分方法找到最优参数值 $\hat{\theta}$：

1. 对 $\hat{\theta_0}$ 求导数。

$\hat{R}(\theta) = \frac{1}{n}\sum^{n}_{i=1} |y_i - \theta|$

$\frac{d}{d\theta} R(\theta) = \frac{d}{d\theta} \left(\frac{1}{n} \sum^{n}_{i=1} |y_i - \theta| \right)$

$= \frac{1}{n} \sum^{n}_{i=1} \frac{d}{d\theta} |y_i - \theta|$

• 这里，我们似乎遇到了一个问题：当参数为 0 时（即 $y_i = \theta$）绝对值的导数是未定义的。现在，我们将忽略这个问题。事实证明，忽略这种情况不会影响我们的最终结果。

• 进行导数运算时，考虑两种情况。当 $\theta$ 小于或等于 $y_i$ 时，项 $y_i - \theta$ 将为正值，绝对值不会产生影响。当 $\theta$ 大于 $y_i$ 时，项 $y_i - \theta$ 将为负值。应用绝对值将其转换为正值，我们可以表示为 $-(y_i - \theta) = \theta - y_i$。

$|y_i - \theta| = \begin{cases} y_i - \theta \quad \text{ 如果 } \theta \le y_i \\ \theta - y_i \quad \text{如果 }\theta > y_i \end{cases}$

• 求导：

$\frac{d}{d\theta} |y_i - \theta| = \begin{cases} \frac{d}{d\theta} (y_i - \theta) = -1 \quad \text{如果 }\theta < y_i \\ \frac{d}{d\theta} (\theta - y_i) = 1 \quad \text{如果 }\theta > y_i \end{cases}$

• 这意味着我们对于 $\theta < y_i$ 和 $\theta > y_i$ 的数据点得到了不同的导数值。我们可以总结为：

$\frac{d}{d\theta} R(\theta) = \frac{1}{n} \sum^{n}_{i=1} \frac{d}{d\theta} |y_i - \theta| \\ = \frac{1}{n} \left[\sum_{\hat{\theta_0} < y_i} (-1) + \sum_{\hat{\theta_0} > y_i} (+1) \right]$

• 换句话说，我们取 $i = 1, 2, ..., n$ 的值的总和：

  • 如果我们的观察值 $y_i$ 大于 我们的预测值 $\hat{\theta_0}$，则为$-1$

  • 如果我们的观察值 $y_i$ 小于 我们的预测值 $\hat{\theta_0}$，则为$+1$

1. 置为 0。$0 = \frac{1}{n}\sum_{\hat{\theta_0} < y_i} (-1) + \frac{1}{n}\sum_{\hat{\theta_0} > y_i} (+1)$

2. 求解 $\hat{\theta_0}$。$0 = -\frac{1}{n}\sum_{\hat{\theta_0} < y_i} (1) + \frac{1}{n}\sum_{\hat{\theta_0} > y_i} (1)$

$\sum_{\hat{\theta_0} < y_i} (1) = \sum_{\hat{\theta_0} > y_i} (1)$

因此，最小化 MAE 的常数模型参数 $\theta = \hat{\theta_0}$ 必须满足：

$\sum_{\hat{\theta_0} < y_i} (1) = \sum_{\hat{\theta_0} > y_i} (1)$

换句话说，大于 $\theta_0$ 的观察数量必须等于小于 $\theta_0$ 的观察数量；方程的左右两侧必须有相等数量的点。这就是中位数的定义，因此我们的最优值是 $\hat{\theta_0} = median(y)$

11.4 总结：损失优化、微积分和临界点

首先，将目标函数定义为平均损失。

• 代入 L1 或 L2 损失。

• 代入模型，使得结果表达为 $\theta$ 的函数。

然后，找到目标函数的最小值：

1. 对 $\theta$ 求导数。

2. 置为 0。

3. 求解 $\hat{\theta}$。

4. （如果我们有多个参数）重复步骤 1-3，使用偏导数。

回想微积分中的临界点：$R(\hat{\theta})$可能是一个最小值、最大值或者鞍点！* 从技术上讲，我们还应该进行二阶导数测试，即，展示 $R''(\hat{\theta}) > 0$。* MSE 具有一个特性——凸性——它保证了 $R(\hat{\theta})$ 是一个全局最小值。* MAE 的凸性证明超出了本课程的范围。

11.5 比较损失函数

我们现在已经尝试了在 MSE 和 MAE 成本函数下拟合模型。这两个结果如何比较？

让我们考虑一个数据集，其中每个条目代表了泡泡茶店每天卖出的饮料数量。我们将拟合一个常数模型来预测明天将卖出的饮料数量。

drinks = np.array([20, 21, 22, 29, 33])
drinks

array([20, 21, 22, 29, 33])

根据我们上面的推导，我们知道 MSE 成本下的最佳模型参数是数据集的均值。在 MAE 成本下，最佳参数是数据集的中位数。

np.mean(drinks), np.median(drinks)

(25.0, 22.0)

如果我们在几个可能的 $\theta$ 值上绘制每个经验风险函数，我们会发现每个 $\hat{\theta}$ 确实对应于最低的错误值：

注意上面的 MSE 是一个平滑函数——它在所有点上都是可微的，这使得用数值方法最小化它变得容易。相比之下，MAE 在每个“拐点”处都不可微。我们将在几周内探讨成本函数的平滑性如何影响我们应用数值优化的能力。

异常值如何影响每个成本函数？想象一下，我们用 1000 替换数据集中的最大值。数据的均值显著增加，而中位数几乎不受影响。

drinks_with_outlier = np.append(drinks, 1033)
display(drinks_with_outlier)
np.mean(drinks_with_outlier), np.median(drinks_with_outlier)

array([  20,   21,   22,   29,   33, 1033])

(193.0, 25.5)

这意味着在 MSE 下，最佳模型参数 $\hat{\theta}$ 受到异常值的影响。在 MAE 下，最佳参数不受异常数据的影响。我们可以通过说 MSE 对异常值敏感，而 MAE 对异常值稳健来概括这一点。

让我们尝试另一个实验。这一次，我们将向数据中添加一个额外的非异常数据点。

drinks_with_additional_observation = np.append(drinks, 35)
drinks_with_additional_observation

array([20, 21, 22, 29, 33, 35])

当我们再次可视化成本函数时，我们发现 MAE 现在在 22 和 29 之间绘制了一条水平线。这意味着模型参数有无数个最佳值：任何值 $\hat{\theta} \in [22, 29]$ 都将最小化 MAE。相比之下，MSE 仍然有一个最佳的 $\hat{\theta}$ 值。换句话说，MSE 有一个唯一的 $\hat{\theta}$ 解；MAE 不能保证有一个唯一的解。

总结我们的例子，

–	MSE（均方损失）	MAE（平均绝对损失）
损失函数	$\hat{R}(\theta) = \frac{1}{n}\sum^{n}_{i=1} (y_i - \theta_0)^2$	$\hat{R}(\theta) = \frac{1}{n}\sum^{n}_{i=1} \|;y_i - \theta_0\|$
最佳 $\hat{\theta_0}$	$\hat{\theta_0} = mean(y) = \bar{y}$	$\hat{\theta_0} = median(y)$
损失曲面
形状	平滑 - 容易使用数值方法最小化（在几周内）	分段 - 在每个“拐点”处，它不可微。更难最小化。
异常值	对异常值敏感（因为它们会显著改变均值）。敏感性还取决于数据集的大小。	对异常值更加稳健。
$\hat{\theta_0}$ 唯一性	唯一 $\hat{\theta_0}$	无数个 $\hat{\theta_0}$

11.6 转换以拟合线性模型

到目前为止，我们已经有了一种有效的方法来拟合模型以预测线性关系。给定一个特征变量和目标，我们可以应用我们的四步过程来找到最佳的模型参数。

上面的关键词是线性。当我们之前计算参数估计时，我们假设$x_i$和$y_i$之间存在大致线性的关系。现实世界中的数据并不总是那么简单，但我们可以对数据进行转换以尝试获得线性关系。

Tukey-Mosteller Bulge Diagram是一个总结两个变量之间关系线性化的变换的有用工具。要确定哪些变换可能合适，追踪数据形成的“凸起”的形状。找到与此凸起匹配的图表象限。该象限的垂直和水平轴上显示的变换可以帮助改善变量之间的拟合。

注意：

• 有多种解决方案。有些比其他的拟合效果更好。

• sqrt 和 log 使值“变小”。

• 提高到幂会使值“变大”。

• 这些变换中的每一个都等同于增加或减少轴的比例。

除了线性之外，还有其他可能的目标，例如使数据看起来更对称。线性允许我们对转换后的数据进行拟合。

让我们重新看一下我们的儒艮示例。长度和年龄如下图所示：

代码

# `corrcoef` computes the correlation coefficient between two variables
# `std` finds the standard deviation
x = dugongs["Length"]
y = dugongs["Age"]
r = np.corrcoef(x, y)[0, 1]
theta_1 = r*np.std(y)/np.std(x)
theta_0 = np.mean(y) - theta_1*np.mean(x)

fig, ax = plt.subplots(1, 2, dpi=200, figsize=(8, 3))
ax[0].scatter(x, y)
ax[0].set_xlabel("Length")
ax[0].set_ylabel("Age")

ax[1].scatter(x, y)
ax[1].plot(x, theta_0 + theta_1*x, "tab:red")
ax[1].set_xlabel("Length")
ax[1].set_ylabel("Age");

在左边的图中，我们看到数据点有轻微的曲线。在右边绘制 SLR 曲线会导致拟合效果不佳。

为了使 SLR 表现良好，我们希望“年龄”和“长度”之间存在粗略的线性趋势。是什么导致原始数据偏离线性关系？注意到“长度”大于 2.6 的数据点相对于其他数据有着不成比例的高“年龄”值。如果我们能够操纵这些数据点使其具有较低的“年龄”值，我们将“移动”这些点向下并减少数据中的曲率。对$y_i$应用对数变换（即取$\log($“年龄”$)$）就可以实现这一点。

关于$\log$的重要说明：在 Data 100（以及大多数高年级 STEM 课程）中，$\log$表示以$e$为底的自然对数。在相关情况下，以 10 为底的对数用$\log_{10}$表示。

代码

z = np.log(y)

r = np.corrcoef(x, z)[0, 1]
theta_1 = r*np.std(z)/np.std(x)
theta_0 = np.mean(z) - theta_1*np.mean(x)

fig, ax = plt.subplots(1, 2, dpi=200, figsize=(8, 3))
ax[0].scatter(x, z)
ax[0].set_xlabel("Length")
ax[0].set_ylabel(r"$\log{(Age)}$")

ax[1].scatter(x, z)
ax[1].plot(x, theta_0 + theta_1*x, "tab:red")
ax[1].set_xlabel("Length")
ax[1].set_ylabel(r"$\log{(Age)}$")

plt.subplots_adjust(wspace=0.3);

我们的 SLR 拟合看起来好多了！我们现在有了一个新的目标变量：SLR 模型现在试图预测“年龄”的对数，而不是未经转换的“年龄”。换句话说，我们应用了变换$z_i = \log{(y_i)}$。注意到得到的模型仍然是参数线性的$\theta = [\theta_0, \theta_1]$。SLR 模型变为：

$\log{\hat{(y_i)}} = \theta_0 + \theta_1 x_i$

$\hat{z}_i = \theta_0 + \theta_1 x_i$

事实证明，这种线性化关系可以帮助我们理解$x_i$和$y_i$之间的基本关系。如果我们重新排列上面的关系，我们会发现：$\log{(y_i)} = \theta_0 + \theta_1 x_i \\ y_i = e^{\theta_0 + \theta_1 x_i} \\ y_i = (e^{\theta_0})e^{\theta_1 x_i} \\ y_i = C e^{k x_i}$

对于一些常数$C$和$k$。

$y_i$是$x_i$的指数函数。对未经转换的变量应用指数拟合可以证实这一发现。

代码

plt.figure(dpi=120, figsize=(4, 3))

plt.scatter(x, y)
plt.plot(x, np.exp(theta_0)*np.exp(theta_1*x), "tab:red")
plt.xlabel("Length")
plt.ylabel("Age");

你可能会想：为什么我们选择特别应用对数变换？为什么不使用其他函数来线性化数据？

实际上，许多其他修改“年龄”和“长度”相对比例的数学运算在这里都可以起作用。

十二、普通最小二乘法

原文：Ordinary Least Squares

译者：飞龙

协议：CC BY-NC-SA 4.0

学习成果

• 定义关于参数向量 $\theta$ 的线性性。

• 了解使用矩阵表示法来表达多元线性回归。

• 解释普通最小二乘法为残差向量的范数的最小化。

• 计算多元线性回归的性能指标。

我们现在已经花了很多讲座时间来探讨如何构建有效的模型 - 我们介绍了 SLR 和常数模型，选择了适合我们建模任务的成本函数，并应用了转换来改进线性拟合。

在所有这些情况下，我们考虑了一个特征的模型 ($\hat{y}_i = \theta_0 + \theta_1 x_i$) 或零个特征的模型 ($\hat{y}_i = \theta_0$)。作为数据科学家，我们通常可以访问包含 许多 特征的数据集。为了建立最佳模型，考虑所有可用的变量作为模型的输入将是有益的，而不仅仅是一个。在今天的讲座中，我们将介绍 多元线性回归 作为将多个特征合并到模型中的框架。我们还将学习如何加速建模过程 - 具体来说，我们将看到线性代数为我们提供了一组强大的工具，用于理解模型性能。

12.1 线性

如果一个表达式是 关于 $\theta$ (一组参数) 是线性组合，那么它是线性的。检查一个表达式是否可以分解为两个项的矩阵乘积 - 一个 $\theta$ 向量，和一个不涉及 $\theta$ 的矩阵/向量 - 是线性的一个很好的指标。

例如，考虑向量 $\theta = [\theta_0, \theta_1, \theta_2]$

1. $\hat{y} = \theta_0 + 2\theta_1 + 3\theta_2$ 在 theta 上是线性的，我们可以将其分解为两个项的矩阵乘积：

$\hat{y} = \begin{bmatrix} 1 \space 2 \space 3 \end{bmatrix} \begin{bmatrix} \theta_0 \\ \theta_1 \\ \theta_2 \end{bmatrix}$

1. $\hat{y} = \theta_0\theta_1 + 2\theta_1^2 + 3log(\theta_2)$ 在 theta 上 不 是线性的，因为 $\theta_1$ 项是平方的，而 $\theta_2$ 项是对数的。我们无法将其分解为两个项的矩阵乘积。

12.2 多元线性回归的术语

在回归的背景下有几个等效的术语。我们在本课程中最常用的是加粗的。

• $x$ 可以被称为

  • 特征

  • 协变量

  • 自变量

  • 解释变量

  • 预测变量

  • 输入

  • 回归器

• $y$ 可以被称为

  • 输出

  • 结果

  • 响应

  • 因变量

• $\hat{y}$ 可以被称为

  • 预测

  • 预测响应

  • 估计值

• $\theta$ 可以被称为

  • 权重

  • 参数

  • 系数

• $\hat{\theta}$ 可以被称为

  • 估计器

  • 最佳参数

• 一个数据点 $(x, y)$ 也被称为一个观测。

12.3 多元线性回归

多元线性回归是简单线性回归的扩展，它将额外的特征添加到模型中。多元线性回归模型的形式为：

$\hat{y} = \theta_0\:+\:\theta_1x_{1}\:+\:\theta_2 x_{2}\:+\:...\:+\:\theta_p x_{p}$

我们对 $y$ 的预测值 $\hat{y}$ 是单个 观测 (特征) $x_i$ 和参数 $\theta_i$ 的线性组合。

我们可以通过查看从 2018-19 NBA 赛季下载的包含每个球员数据的数据集来进一步探讨这个想法，数据来自 Kaggle。

代码

import pandas as pd
nba = pd.read_csv('data/nba18-19.csv', index_col=0)
nba.index.name = None # Drops name of index (players are ordered by rank)

nba.head(5)

	Player	Pos	Age	Tm	G	GS	MP	FG	FGA	FG%	...	FT%	ORB	DRB	TRB	AST	STL	BLK	TOV	PF	PTS
1	Álex Abrines\abrinal01	SG	25	OKC	31	2	19.0	1.8	5.1	0.357	...	0.923	0.2	1.4	1.5	0.6	0.5	0.2	0.5	1.7	5.3
2	Quincy Acy\acyqu01	PF	28	PHO	10	0	12.3	0.4	1.8	0.222	...	0.700	0.3	2.2	2.5	0.8	0.1	0.4	0.4	2.4	1.7
3	Jaylen Adams\adamsja01	PG	22	ATL	34	1	12.6	1.1	3.2	0.345	...	0.778	0.3	1.4	1.8	1.9	0.4	0.1	0.8	1.3	3.2
4	Steven Adams\adamsst01	C	25	OKC	80	80	33.4	6.0	10.1	0.595	...	0.500	4.9	4.6	9.5	1.6	1.5	1.0	1.7	2.6	13.9
5	Bam Adebayo\adebaba01	C	21	MIA	82	28	23.3	3.4	5.9	0.576	...	0.735	2.0	5.3	7.3	2.2	0.9	0.8	1.5	2.5	8.9

5 行×29 列

假设我们有兴趣预测一名运动员本赛季在篮球比赛中将得分的数量（PTS）。

假设我们想要通过使用球员的一些特征或特征来拟合一个线性模型。具体来说，我们将专注于投篮命中、助攻和三分球出手。

• FG，每场比赛的（2 分）投篮命中数

• AST，每场比赛的平均助攻数

• 3PA，每场比赛尝试的三分球数

代码

nba[['FG', 'AST', '3PA', 'PTS']].head()

	FG	AST	3PA	PTS
1	1.8	0.6	4.1	5.3
2	0.4	0.8	1.5	1.7
3	1.1	1.9	2.2	3.2
4	6.0	1.6	0.0	13.9
5	3.4	2.2	0.2	8.9

因为我们现在处理的是许多参数值，我们已经将它们全部收集到了一个维度为$(p+1) \times 1$的参数向量中，以保持整洁。记住$p$代表我们拥有的特征数量（在这种情况下是 3）。

$\theta = \begin{bmatrix} \theta_{0} \\ \theta_{1} \\ \vdots \\ \theta_{p} \end{bmatrix}$

我们在这里使用两个向量：一个表示观察数据的行向量，另一个包含模型参数的列向量。多元线性回归模型等同于观察向量和参数向量的点（标量）积。

$[1,\:x_{1},\:x_{2},\:x_{3},\:...,\:x_{p}] \theta = [1,\:x_{1},\:x_{2},\:x_{3},\:...,\:x_{p}] \begin{bmatrix} \theta_{0} \\ \theta_{1} \\ \vdots \\ \theta_{p} \end{bmatrix} = \theta_0\:+\:\theta_1x_{1}\:+\:\theta_2 x_{2}\:+\:...\:+\:\theta_p x_{p}$

请注意，我们已经在观察向量中插入了 1 作为第一个值。当计算点积时，这个 1 将与$\theta_0$相乘，得到回归模型的截距。我们称这个 1 条目为截距或偏差项。

鉴于我们这里有三个特征，我们可以将这个模型表示为：$\hat{y} = \theta_0\:+\:\theta_1x_{1}\:+\:\theta_2 x_{2}\:+\:\theta_3 x_{3}$

我们的特征由$x_1$（FG）、$x_2$（AST）和$x_3$（3PA）表示，每个特征都有对应的参数$\theta_1$、$\theta_2$和$\theta_3$。

在统计学中，这个模型+损失被称为普通最小二乘法（OLS）。OLS 的解是参数$\hat{\theta}$的最小损失，也称为最小二乘估计。

12.4 线性代数方法

我们现在知道如何从多个观察特征生成单个预测。数据科学家通常会进行大规模工作 - 也就是说，他们希望构建可以一次产生多个预测的模型。我们上面介绍的向量表示法为我们提供了如何加速多元线性回归的线索。我们想要使用线性代数的工具。

让我们考虑如何应用上面所做的事情。为了适应我们正在考虑多个特征变量的事实，我们将稍微调整我们的符号。现在，每个观察可以被认为是一个行向量，其中每个特征都有一个条目。

要从数据中的第一个观测中进行预测，我们取参数向量和第一个观测向量的点积。要从第二个观测中进行预测，我们将重复这个过程，找到参数向量和第二个观测向量的点积。如果我们想要找到数据集中每个观测的模型预测，我们将对数据中的所有$n$个观测重复这个过程。

$\hat{y}_1 = \theta_0 + \theta_1 x_{11} + \theta_2 x_{12} + ... + \theta_p x_{1p} = [1,\:x_{11},\:x_{12},\:x_{13},\:...,\:x_{1p}] \theta$

$\hat{y}_2 = \theta_0 + \theta_1 x_{21} + \theta_2 x_{22} + ... + \theta_p x_{2p} = [1,\:x_{21},\:x_{22},\:x_{23},\:...,\:x_{2p}] \theta$

$\vdots$

$\hat{y}_n = \theta_0 + \theta_1 x_{n1} + \theta_2 x_{n2} + ... + \theta_p x_{np} = [1,\:x_{n1},\:x_{n2},\:x_{n3},\:...,\:x_{np}] \theta$

我们的观测数据由$n$个行向量表示，每个向量的维度为$(p+1)$。我们可以将它们全部收集到一个称为$\mathbb{X}$的单个矩阵中。

矩阵$\mathbb{X}$被称为设计矩阵。它包含了我们$p$个特征的所有观测数据，其中每一行对应一个观测，每一列对应一个特征。它通常（但并非总是）包含一个额外的全为 1 的列来表示截距或偏置列。

回顾设计矩阵中发生的情况：每一行代表一个单独的观测。例如，数据 100 中的一个学生。每一列代表一个特征。例如，数据 100 中学生的年龄。这个约定使我们能够轻松地将我们在数据框中的先前工作转移到这种新的线性代数视角。

多元线性回归模型可以用矩阵的术语重新表述：$\Large \mathbb{\hat{Y}} = \mathbb{X} \theta$

在这里，$\mathbb{\hat{Y}}$是具有$n$个元素的预测向量（$\mathbb{\hat{Y}} \in \mathbb{R}^{n}$）；它包含模型对每个$n$个输入观测的预测。 $\mathbb{X}$是维度为$\mathbb{X} \in \mathbb{R}^(n \times (p + 1))$的设计矩阵，$\theta$是维度为$\theta \in \mathbb{R}^{(p + 1)}$的参数向量。

作为一个复习，让我们也回顾一下点积（或内积）。这是一个向量运算，它：

• 只能在两个相同长度的向量上进行

• 对应向量的乘积求和

• 返回一个单一的数字

虽然这不在范围内，但请注意我们也可以几何地解释点积：

• 它是三个因素的乘积：两个向量的大小和它们之间的角度的余弦：$\vec{u} \cdot \vec{v} = ||\vec{u}|| \cdot ||\vec{v}|| \cdot {cos \theta}$

12.5 均方误差

现在我们有了一个新的理解模型的方法，以向量和矩阵为基础。为了配合这个新的约定，我们应该更新我们对风险函数和模型拟合的理解。

回想一下我们对 MSE 的定义：$R(\theta) = \frac{1}{n} \sum_{i=1}^n (y_i - \hat{y}_i)^2$

在本质上，MSE 是一种距离的度量 - 它指示了预测值与真实值之间的“距离”平均有多远。

在处理向量时，这种“距离”或向量的大小/长度的概念由范数表示。更确切地说，向量$\vec{a}$和$\vec{b}$之间的距离可以表示为：$||\vec{a} - \vec{b}||_2 = \sqrt{(a_1 - b_1)^2 + (a_2 - b_2)^2 + \ldots + (a_n - b_n)^2} = \sqrt{\sum_{i=1}^n (a_i - b_i)^2}$

双竖线是范数的数学表示。下标 2 表示我们正在计算 L2 范数，或平方范数。

我们需要了解 Data 100 的两种范数是 L1 和 L2 范数（听起来熟悉吗？）。在这篇笔记中，我们将专注于 L2 范数。我们将在未来的讲座中深入探讨 L1 范数。

对于 n 维向量$\vec{x} = \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{bmatrix}$，L2 向量范数是

$||\vec{x}||_2 = \sqrt{(x_1)^2 + (x_2)^2 + \ldots + (x_n)^2} = \sqrt{\sum_{i=1}^n (x_i)^2}$

L2 向量范数是$n$维中勾股定理的推广。因此，它可以用作矢量的长度的度量，甚至可以用作两个矢量之间的距离的度量。

我们可以将 MSE 表示为平方 L2 范数，如果我们用预测向量$\hat{\mathbb{Y}}$和真实目标向量$\mathbb{Y}$来重新表达它：

$R(\theta) = \frac{1}{n} \sum_{i=1}^n (y_i - \hat{y}_i)^2 = \frac{1}{n} ||\mathbb{Y} - \hat{\mathbb{Y}}||_2^2$

这里，范数双条之外的上标 2 表示我们正在平方范数。如果我们插入我们的线性模型$\hat{\mathbb{Y}} = \mathbb{X} \theta$，我们会发现 MSE 成本函数的向量表示：

$R(\theta) = \frac{1}{n} ||\mathbb{Y} - \mathbb{X} \theta||_2^2$

在线性代数的视角下，我们的新任务是拟合最佳参数向量$\theta$，使得成本函数最小化。等价地，我们希望最小化范数$||\mathbb{Y} - \mathbb{X} \theta||_2 = ||\mathbb{Y} - \hat{\mathbb{Y}}||_2.$

我们可以用两种方式重新表述这个目标：

• 最小化真实值向量$\mathbb{Y}$和预测值向量$\mathbb{\hat{Y}}$之间的距离

• 最小化残差向量的长度，定义为：$e = \mathbb{Y} - \mathbb{\hat{Y}} = \begin{bmatrix} y_1 - \hat{y}_1 \\ y_2 - \hat{y}_2 \\ \vdots \\ y_n - \hat{y}_n \end{bmatrix}$

12.6 几何视角

为了得出最佳参数向量以实现这一目标，我们可以利用我们建模设置的几何特性。

到目前为止，我们大多把我们的模型看作是观测值和参数向量水平堆叠的标量积。我们也可以将$\hat{\mathbb{Y}}$看作是特征向量的线性组合，由参数缩放。我们使用符号$\mathbb{X}_{:, i}$来表示设计矩阵的第$i$列。您可以将其视为在调用.iloc和.loc时使用的相同约定。“:”表示我们正在取第$i$列中的所有条目。

$\hat{\mathbb{Y}} = \theta_0 \begin{bmatrix} 1 \\ 1 \\ \vdots \\ 1 \end{bmatrix} + \theta_1 \begin{bmatrix} x_{11} \\ x_{21} \\ \vdots \\ x_{n1} \end{bmatrix} + \ldots + \theta_p \begin{bmatrix} x_{1p} \\ x_{2p} \\ \vdots \\ x_{np} \end{bmatrix} = \theta_0 \mathbb{X}_{:,\:1} + \theta_1 \mathbb{X}_{:,\:2} + \ldots + \theta_p \mathbb{X}_{:,\:p+1}$

这种新方法很有用，因为它使我们能够利用线性组合的性质。

回想一下矩阵$\mathbb{X}$的范围或列空间（表示为$span(\mathbb{X})$）是矩阵列的所有可能线性组合的集合。换句话说，范围代表着可能通过添加和缩放矩阵列的某些组合到达的空间中的每一点。另外，如果$\mathbb{X}$的每一列的长度为$n$，$span(\mathbb{X})$是$\mathbb{R}^{n}$的子空间。

因为预测向量$\hat{\mathbb{Y}} = \mathbb{X} \theta$是$\mathbb{X}$的列的线性组合，我们知道预测包含在$\mathbb{X}$的范围内。也就是说，我们知道$\mathbb{\hat{Y}} \in \text{Span}(\mathbb{X})$。

下面的图是对$\text{Span}(\mathbb{X})$的简化视图，假设$\mathbb{X}$的每一列都有长度$n$。注意\mathbb{X的列定义了$\mathbb{R}^n$的子空间，子空间中的每个点都可以通过$\mathbb{X}$的列的线性组合到达。预测向量$\mathbb{\hat{Y}}$位于这个子空间的某个位置。

检查这个图，我们发现了一个问题。真实值向量$\mathbb{Y}$理论上可以位于$\mathbb{R}^n$空间中的任何位置——它的确切位置取决于我们在现实世界中收集的数据。然而，我们的多元线性回归模型只能在$\mathbb{X}$张成的$\mathbb{R}^n$空间的子空间中进行预测。记住我们在前一节建立的模型拟合目标：我们希望生成预测，使得真实值向量$\mathbb{Y}$和预测值向量$\mathbb{\hat{Y}}$之间的距离最小化。这意味着我们希望$\mathbb{\hat{Y}}$是$\text{Span}(\mathbb{X})$中离$\mathbb{Y}$最近的向量。

另一种重新表述这个目标的方式是，我们希望最小化残差向量$e$的长度，即其$L_2$范数。

在$\text{Span}(\mathbb{X})$中距离$\mathbb{Y}$最近的向量始终是$\mathbb{Y}$在$\text{Span}(\mathbb{X})$上的正交投影。因此，我们应该选择参数向量$\theta$，使得残差向量与$\text{Span}(\mathbb{X})$中的任何向量正交。你可以将这个想象成从$\mathbb{Y}$到$\mathbb{X}$的跨度上垂直投影线创建的向量。

这如何帮助我们确定最佳参数向量$\hat{\theta}$？回想一下，如果两个向量$a$和$b$正交，它们的点积为零：${a}^{T}b = 0$。如果向量$v$正交于矩阵$M$的张成空间，当且仅当$v$正交于$M$中的每一列。综合起来，向量$v$对于$\text{Span}(M)$正交，如果：

$M^Tv = \vec{0}$

请注意，$\vec{0}$代表零向量，一个全为 0 的$d$长度向量。

记住我们的目标是找到$\hat{\theta}$，使得我们最小化目标函数$R(\theta)$。等价地，这就是使得残差向量$e = \mathbb{Y} - \mathbb{X} \theta$与$\text{Span}(\mathbb{X})$正交的$\hat{\theta}$。

观察$\mathbb{Y} - \mathbb{X}\hat{\theta}$与$span(\mathbb{X})$正交的定义（0 是$\vec{0}$向量），我们可以写成：$\mathbb{X}^T (\mathbb{Y} - \mathbb{X}\hat{\theta}) = \vec{0}$

然后我们重新排列项：$\mathbb{X}^T \mathbb{Y} - \mathbb{X}^T \mathbb{X} \hat{\theta} = \vec{0}$

最后，我们得到了正规方程：$\mathbb{X}^T \mathbb{X} \hat{\theta} = \mathbb{X}^T \mathbb{Y}$

任何最小化数据集上均方误差的向量$\theta$必须满足这个方程。

如果$\mathbb{X}^T \mathbb{X}$是可逆的，我们可以得出结论：$\hat{\theta} = (\mathbb{X}^T \mathbb{X})^{-1} \mathbb{X}^T \mathbb{Y}$

这被称为$\theta$的最小二乘估计：它是使平方损失最小化的$\theta$的值。

请注意，最小二乘估计是在假设$\mathbb{X}^T \mathbb{X}$是可逆的条件下推导出来的。当$\mathbb{X}^T \mathbb{X}$是满列秩时，这个条件成立，而这又发生在$\mathbb{X}$是满列秩时。我们将在实验和作业中探讨这个事实的后果。

12.7 评估模型性能

我们对多元线性回归的几何视图已经有了很大的进展！我们已经确定了最小化多个特征模型中的均方误差的参数值的最佳集合。

现在，我们想要了解我们的拟合模型的表现如何。模型性能的一个度量是均方根误差，即 RMSE。RMSE 只是 MSE 的平方根。取平方根将值转换回$y_i$的原始、非平方单位，这对于理解模型的性能很有用。较低的 RMSE 表示更“准确”的预测-在整个数据集中有更低的平均损失。

$\text{RMSE} = \sqrt{\frac{1}{n} \sum_{i=1}^n (y_i - \hat{y}_i)^2}$

在处理 SLR 时，我们生成了残差与单个特征的图表，以了解残差的行为。在多元线性回归中使用多个特征时，考虑在残差图中只有一个特征不再有意义。相反，多元线性回归通过制作残差与预测值的图表来进行评估。与 SLR 一样，如果多元线性模型的残差图没有模式，则表现良好。

对于 SLR，我们使用相关系数来捕捉目标变量和单个特征变量之间的关联。在多元线性模型设置中，我们将需要一个性能度量，可以同时考虑多个特征。多元$R^2$，也称为决定系数，是我们的拟合值（预测）$\hat{y}_i$的方差比例到真实值$y_i$。它的范围从 0 到 1，实际上是模型解释观察中方差的比例。

$R^2 = \frac{\text{variance of } \hat{y}_i}{\text{variance of } y_i} = \frac{\sigma^2_{\hat{y}}}{\sigma^2_y}$

请注意，对于具有截距项的 OLS，例如$\hat{y} = \theta_0 + \theta_1x_1 + \theta_2x_2 + \cdots + \theta_px_p$，$\mathbb{R}^2$等于$y$和$\hat{y}$之间的相关性的平方。另一方面，对于 SLR，$\mathbb{R}^2$等于$r^2$，即$x$和$y$之间的相关性。这两个属性的证明超出了本课程的范围。

此外，随着我们添加更多的特征，我们的拟合值倾向于越来越接近我们的实际值。因此，$\mathbb{R}^2$增加。

然而，增加更多的特征并不总是意味着我们的模型更好！我们将在课程后面看到原因。

12.8 OLS 属性

1. 使用最优参数向量时，我们的残差$e = \mathbb{Y} - \hat{\mathbb{Y}}$与$span(\mathbb{X})$正交。

$\mathbb{X}^Te = 0$

证明：

• 最优参数向量$\hat{\theta}$解决了正规方程$\implies \hat{\theta} = \mathbb{X}^T\mathbb{X}^{-1}\mathbb{X}^T\mathbb{Y}$

$\mathbb{X}^Te = \mathbb{X}^T (\mathbb{Y} - \mathbb{\hat{Y}})$

$\mathbb{X}^T (\mathbb{Y} - \mathbb{X}\hat{\theta}) = \mathbb{X}^T\mathbb{Y} - \mathbb{X}^T\mathbb{X}\hat{\theta}$

• 任何矩阵与其逆矩阵相乘都是单位矩阵$\mathbb{I}$

$\mathbb{X}^T\mathbb{Y} - (\mathbb{X}^T\mathbb{X})(\mathbb{X}^T\mathbb{X})^{-1}\mathbb{X}^T\mathbb{Y} = \mathbb{X}^T\mathbb{Y} - \mathbb{X}^T\mathbb{Y} = 0$

2. 对于所有具有截距项的线性模型，残差的总和为零。

$\sum_i^n e_i = 0$

证明：

• 对于所有具有截距项的线性模型，预测的$y$值的平均值等于真实$y$值的平均值。$\bar{y} = \bar{\hat{y}}$

• 将残差总和重写为两个单独的总和，$\sum_i^n e_i = \sum_i^n y_i - \sum_i^n\hat{y}_i$

• 每个相应的和是平均和的倍数。$\sum_i^n e_i = n\bar{y} - n\bar{y} = n(\bar{y} - \bar{y}) = 0$

3. 最小二乘估计$\hat{\theta}$是唯一的，当且仅当$\mathbb{X}$是满列秩的。

证明：

• 我们知道正规方程$\mathbb{X}^T\mathbb{X}\hat{\theta} = \mathbb{Y}$的解是满足先前相等的最小二乘估计。

• $\hat{\theta}$ 有一个唯一的解 $\iff$ 方阵 $\mathbb{X}^T\mathbb{X}$ 是可逆的。

  • 方阵的列秩是它包含的线性独立列的数量。

  • 一个 $n$ x $n$ 的方阵被认为是完整的列秩当且仅当它的所有列都是线性独立的。也就是说，它的秩等于 $n$。

  • $\mathbb{X}^T\mathbb{X}$ 的形状是 $(p + 1) \times (p + 1)$，因此最大秩为 $p + 1$。

• $rank(\mathbb{X}^T\mathbb{X})$ = $rank(\mathbb{X})$（证明超出范围）。

• 因此，$\mathbb{X}^T\mathbb{X}$ 的秩为 $p + 1$ $\iff$ $\mathbb{X}$ 的秩为 $p + 1$ $\iff \mathbb{X}$ 是完整的列秩。

总结：

	模型	估计	唯一？
常数模型 + MSE	$\hat{y} = \theta_0$	$\hat{\theta_0} = mean(y) = \bar{y}$	是。任何一组值都有唯一的均值。
常数模型 + MAE	$\hat{y} = \theta_0$	$\hat{\theta_0} = median(y)$	是，如果是奇数。否，如果是偶数。返回中间 2 个值的平均值。
简单线性回归 + MSE	$\hat{y} = \theta_0 + \theta_1x$	$\hat{\theta_0} = \bar{y} - \hat{\theta_1}\hat{x}$ $\hat{\theta_1} = r\frac{\sigma_y}{\sigma_x}$	是。任何一组非常数*值都有唯一的均值、标准差和相关系数。
OLS（线性模型 + MSE）	$\mathbb{\hat{Y}} = \mathbb{X}\mathbb{\theta}$	$\hat{\theta} = \mathbb{X}^T\mathbb{X}^{-1}\mathbb{X}^T\mathbb{Y}$	是，如果 $\mathbb{X}$ 是完整的列秩（所有列线性独立，数据点的数量 >>> 特征的数量）。

十三、梯度下降

原文：Gradient Descent

译者：飞龙

协议：CC BY-NC-SA 4.0

学习成果

• 优化复杂模型

• 识别直接微积分或几何论证无法帮助解决损失函数的情况

• 应用梯度下降进行数值优化

到目前为止，我们已经非常熟悉选择模型和相应损失函数的过程，并通过选择最小化损失函数的$\theta$的值来优化参数。到目前为止，我们已经通过以下两种方法优化了$\theta$：1. 使用微积分对损失函数关于$\theta$的导数进行求导，将其置为 0，并解出$\theta$。2. 使用正交性的几何论证来推导 OLS 解$\hat{\theta} = (\mathbb{X}^T \mathbb{X})^{-1}\mathbb{X}^T \mathbb{Y}$。

然而，需要注意的一点是，我们上面使用的技术只有在我们做出一些重大假设时才能应用。对于微积分方法，我们假设损失函数在所有点上都是可微的，并且代数是可管理的；对于几何方法，OLS仅适用于使用 MSE 损失的线性模型。当我们有更复杂的模型和不同的更复杂的损失函数时会发生什么？到目前为止我们学到的技术将不起作用，所以我们需要一种新的优化技术：梯度下降。

重要思想：使用算法而不是求解精确答案

13.1 最小化 1D 函数

让我们考虑一个任意的函数。我们的目标是找到最小化这个函数的$x$的值。

代码

import pandas as pd
import seaborn as sns
import plotly.express as px
import matplotlib.pyplot as plt
import numpy as np
pd.options.mode.chained_assignment = None  # default='warn'

def arbitrary(x):
 return (x**4 - 15*x**3 + 80*x**2 - 180*x + 144)/10

13.1.1 天真的方法：猜测和检查

以上，我们看到最小值大约在 5.3 左右。让我们看看我们是否能够弄清楚如何从头开始算法地找到确切的最小值。一种非常慢（而且糟糕）的方法是手动猜测和检查。

arbitrary(6)

0.0

一个稍微好一点（但仍然很慢）的方法是使用蛮力来尝试一堆 x 值，并返回产生最低损失的值。

def simple_minimize(f, xs):
 # Takes in a function f and a set of values xs. 
 # Calculates the value of the function f at all values x in xs
 # Takes the minimum value of f(x) and returns the corresponding value x 
 y = [f(x) for x in xs] 
 return xs[np.argmin(y)]

guesses = [5.3, 5.31, 5.32, 5.33, 5.34, 5.35]
simple_minimize(arbitrary, guesses)

5.33

这个过程本质上与我们以前制作图形图表的过程相同，只是我们只看了 20 个选定的点。

代码

xs = np.linspace(1, 7, 200)
sparse_xs = np.linspace(1, 7, 5)

ys = arbitrary(xs)
sparse_ys = arbitrary(sparse_xs)

fig = px.line(x = xs, y = arbitrary(xs))
fig.add_scatter(x = sparse_xs, y = arbitrary(sparse_xs), mode = "markers")
fig.update_layout(showlegend= False)
fig.update_layout(autosize=False, width=800, height=600)
fig.show()

这种基本方法存在三个主要缺陷：1. 如果最小值在我们猜测的范围之外，答案将完全错误。2. 即使我们的猜测范围是正确的，如果猜测太粗糙，我们的答案将不准确。3. 考虑到可能庞大的无用猜测数量，这是绝对计算效率低下的。

13.1.2 Scipy.optimize.minimize

最小化这个数学函数的一种方法是使用scipy.optimize.minimize函数。它接受一个函数和一个起始猜测，并尝试找到最小值。

from scipy.optimize import minimize

# takes a function f and a starting point x0 and returns a readout 
# with the optimal input value of x which minimizes f
minimize(arbitrary, x0 = 3.5)

 message: Optimization terminated successfully.
  success: True
   status: 0
      fun: -0.13827491292966557
        x: [ 2.393e+00]
      nit: 3
      jac: [ 6.486e-06]
 hess_inv: [[ 7.385e-01]]
     nfev: 20
     njev: 10

scipy.optimize.minimize很棒。它也可能看起来有点神奇。你怎么能写一个函数来找到任何数学函数的最小值呢？有许多方法可以做到这一点，我们将在今天的讲座中探讨这些方法，最终到达scipy.optimize.minimize使用的梯度下降的重要思想。

事实证明，在幕后，LinearRegression模型的fit方法使用了梯度下降。梯度下降也是许多机器学习模型的工作原理，甚至包括先进的神经网络模型。

在 Data 100 中，梯度下降过程通常对我们来说是不可见的，隐藏在抽象层下。然而，作为优秀的数据科学家，我们知道优化函数利用的基本原理是很重要的。

13.2 深入研究梯度下降

在这个域中查看函数，很明显函数的最小值出现在$\theta = 5.3$附近。假设一下，我们看不到成本函数的完整视图。我们如何猜测最小化函数的$\theta$值？

原来，函数的一阶导数可以给我们一些线索。在下面的图中，线表示每个$\theta$值的导数值。导数为负值时为红色，为正值时为绿色。

假设我们对最小化$\theta$的值进行了猜测。记住我们从左到右读取图，并假设我们的起始$\theta$值在最佳$\hat{\theta}$的左边。如果猜测“低估”了真正的最小值 - 我们对最小化函数的$\hat{\theta}$的猜测低于真正的$\hat{\theta}$值 - 导数将是负数。这意味着如果我们增加$\theta$（向右移动），那么我们可以进一步减少我们的损失函数。如果这个猜测“高估”了真正的最小值，导数将是正数，暗示着相反。

我们可以利用这种模式来帮助制定我们对最佳$\hat{\theta}$的下一个猜测。考虑一种情况，我们通过猜测一个太低的值而低估了$\theta$。我们希望我们的下一个猜测的值比上一个猜测的值更大 - 也就是说，我们希望将我们的猜测向右移动。你可以把这看作是沿着斜坡“下坡”到函数的最小值。

如果我们通过猜测一个太高的值而高估了$\hat{\theta}$，我们希望我们的下一个猜测的值更低 - 我们希望将$\hat{\theta}$的猜测向左移动。

换句话说，每个点的函数的导数告诉我们我们下一个猜测的方向。负斜率意味着我们想向右走，或者向正方向移动。正斜率意味着我们想向左走，或者向负方向移动。

13.2.1 算法尝试 1

有了这个知识，让我们试着看看我们是否可以使用导数来优化函数。

我们首先对$x$的最小值进行一些猜测。然后，我们查看该$x$值的函数的导数，并向相反方向下坡。我们可以将我们的新规则表示为一个递归关系：

$x^{(t+1)} = x^{(t)} - \frac{d}{dx} f(x^{(t)})$

将这个陈述翻译成英文：我们通过取上一次的猜测($x^{(t)}$)并减去该点的函数的导数($\frac{d}{dx} f(x^{(t)})$)来获得我们下一次的猜测，即在时间步长$t+1$($x^{(t+1)}$)的最小值$x$。

下面显示了一些步骤，其中旧步骤显示为透明点，下一个步骤是填充为绿色的点。

看起来不错！但是我们确实有一个问题 - 一旦我们接近函数的最小值，我们的猜测就会在最小值附近“反弹”而永远无法到达它。

换句话说，每次更新我们的猜测时，我们走得太远了。我们可以通过减小每个步骤的大小来解决这个问题。

13.2.2 算法尝试 2

让我们更新我们的算法以使用学习率（有时也称为步长），它控制我们每次更新的距离。我们用$\alpha$表示学习率。

$x^{(t+1)} = x^{(t)} - \alpha \frac{d}{dx} f(x^{(t)})$

小的$\alpha$意味着我们会采取小步；大的$\alpha$意味着我们会采取大步。

更新我们的函数使用$\alpha=0.3$，我们的算法成功地收敛（在最小值上解决并停止显着更新，或者完全停止）。

13.2.3 凸性

在我们上面的分析中，我们把注意力集中在损失函数的全局最小值上。你可能会想：左边的局部最小值呢？

如果我们选择了不同的$\theta$起始猜测，或者不同的学习率$\alpha$值，我们的算法可能会“卡住”，收敛到局部最小值，而不是真正的最优损失值。

如果损失函数是凸的，梯度下降保证会收敛并找到目标函数的全局最小值。形式上，如果一个函数$f$满足：$tf(a) + (1-t)f(b) \geq f(ta + (1-t)b)$ 对于所有$f$的定义域中的$a, b$和$t \in [0, 1]$。

换句话说：如果你在曲线上的任意两点之间画一条线，曲线上的所有值必须在或下于该线。重要的是，凸函数的任何局部最小值也是它的全局最小值。

总之，非凸损失函数可能会导致优化问题。这意味着我们选择的损失函数是建模过程中的关键因素。事实证明，MSE 是凸的，这也是它成为如此受欢迎的损失函数的主要原因。

13.3 一维梯度下降

术语澄清：在过去的讲座中，我们使用“损失”来指代单个数据点上发生的错误。在应用中，我们通常更关心所有数据点的平均误差。未来，我们将把“模型的损失”理解为数据集上的模型平均误差。这有时也被称为经验风险、成本函数或目标函数。$L(\theta) = R(\theta) = \frac{1}{n} \sum_{i=1}^{n} L(y, \hat{y})$

在上面的讨论中，我们使用了一些任意函数$f$。作为数据科学家，我们几乎总是在优化模型的情况下使用梯度下降 - 具体来说，我们想要应用梯度下降来找到损失函数的最小值。在建模的情况下，我们的目标是通过选择最小化的模型参数来最小化损失函数。

回顾一下我们过去几节讲座的建模工作流程：* 定义具有一些参数$\theta_i$的模型 * 选择一个损失函数 * 选择使数据上的损失函数最小化的$\theta_i$的值

梯度下降是完成最后任务的强大技术。通过应用梯度下降算法，我们可以选择参数$\theta_i$的值，这将导致模型在训练数据上损失最小。

在建模上下文中使用梯度下降时：* 我们对最小化的$\theta_i$进行猜测 * 我们计算损失函数$L$的导数

我们可以通过用$\theta$替换$x$和用$L$替换$f$来“翻译”我们之前的梯度下降规则：

$\theta^{(t+1)} = \theta^{(t)} - \alpha \frac{d}{d\theta} L(\theta^{(t)})$

13.3.1 在tips数据集上的梯度下降

为了看到这一点，让我们考虑一个没有偏移的线性模型的情况。我们想要预测小费（y）给定一顿饭的价格（x）。为了做到这一点，我们

• 选择一个模型：$\hat{y} = \theta_1 x$，

• 选择一个损失函数：$L(\theta) = MSE(\theta) = \frac{1}{n} \sum_{i=1}^n (y_i - \theta_1x_i)^2$。

让我们应用之前的gradient_descent函数来优化我们在tips数据集上的模型。我们将尝试选择最佳参数$\theta_i$来预测total_bill$x$的tip$y$。

df = sns.load_dataset("tips")
df.head()

	total_bill	tip	sex	smoker	day	time	size
0	16.99	1.01	Female	No	Sun	Dinner	2
1	10.34	1.66	Male	No	Sun	Dinner	3
2	21.01	3.50	Male	No	Sun	Dinner	3
3	23.68	3.31	Male	No	Sun	Dinner	2
4	24.59	3.61	Female	No	Sun	Dinner	4

我们可以可视化我们的数据集上 MSE 的值，对于不同可能的$\theta_1$的选择。为了优化我们的模型，我们希望选择导致最低 MSE 的$\theta_1$的值。

代码

import plotly.graph_objects as go

def derivative_arbitrary(x):
 return (4*x**3 - 45*x**2 + 160*x - 180)/10

fig = go.Figure()
roots = np.array([2.3927, 3.5309, 5.3263])

fig.add_trace(go.Scatter(x = xs, y = arbitrary(xs), 
 mode = "lines", name = "f"))
fig.add_trace(go.Scatter(x = xs, y = derivative_arbitrary(xs), 
 mode = "lines", name = "df", line = {"dash": "dash"}))
fig.add_trace(go.Scatter(x = np.array(roots), y = 0*roots, 
 mode = "markers", name = "df = zero", marker_size = 12))
fig.update_layout(font_size = 20, yaxis_range=[-1, 3])
fig.update_layout(autosize=False, width=800, height=600)
fig.show()

要应用梯度下降，我们需要计算损失函数对我们的参数$\theta_1$的导数。

• 给定我们的损失函数，$L(\theta) = MSE(\theta) = \frac{1}{n} \sum_{i=1}^n (y_i - \theta_1x_i)^2$

• 我们对$\theta_1$进行导数$\frac{\partial}{\partial \theta_{1}} L(\theta_1^{(t)}) = \frac{-2}{n} \sum_{i=1}^n (y_i - \theta_1^{(t)} x_i) x_i$

• 这导致梯度下降更新规则$\theta_1^{(t+1)} = \theta_1^{(t)} - \alpha \frac{d}{d\theta}L(\theta_1^{(t)})$

对于一些学习率$\alpha$。

在代码中实现这一点，我们可以可视化tips数据上的 MSE 损失。MSE 是凸的，因此有一个全局最小值。

代码

def gradient_descent(df, initial_guess, alpha, n):
 """Performs n steps of gradient descent on df using learning rate alpha starting
 from initial_guess. Returns a numpy array of all guesses over time."""
 guesses = [initial_guess]
 current_guess = initial_guess
 while len(guesses) < n:
 current_guess = current_guess - alpha * df(current_guess)
 guesses.append(current_guess)

 return np.array(guesses)

def mse_single_arg(theta_1):
 """Returns the MSE on our data for the given theta1"""
 x = df["total_bill"]
 y_obs = df["tip"]
 y_hat = theta_1 * x
 return np.mean((y_hat - y_obs) ** 2)

def mse_loss_derivative_single_arg(theta_1):
 """Returns the derivative of the MSE on our data for the given theta1"""
 x = df["total_bill"]
 y_obs = df["tip"]
 y_hat = theta_1 * x

 return np.mean(2 * (y_hat - y_obs) * x)

loss_df = pd.DataFrame({"theta_1":np.linspace(-1.5, 1), "MSE":[mse_single_arg(theta_1) for theta_1 in np.linspace(-1.5, 1)]})

trajectory = gradient_descent(mse_loss_derivative_single_arg, -0.5, 0.0001, 100)

plt.plot(loss_df["theta_1"], loss_df["MSE"])
plt.scatter(trajectory, [mse_single_arg(guess) for guess in trajectory], c="white", edgecolor="firebrick")
plt.scatter(trajectory[-1], mse_single_arg(trajectory[-1]), c="firebrick")
plt.xlabel(r"$\theta_1$")
plt.ylabel(r"$L(\theta_1)$");

print(f"Final guess for theta_1: {trajectory[-1]}")

Final guess for theta_1: 0.14369554654231262

13.4 多维模型的梯度下降

我们上面使用的函数是一维的 - 我们只是在最小化函数相对于单个参数$\theta$。然而，模型通常具有需要优化的多个参数的成本函数。例如，简单线性回归有 2 个参数：$\hat{y} + \theta_0 + \theta_1x$，多元线性回归有$p+1$个参数：$\mathbb{Y} = \theta_0 + \theta_1 \Bbb{X}_{:,1} + \theta_2 \Bbb{X}_{:,2} + \cdots + \theta_p \Bbb{X}_{:,p}$

我们需要扩展梯度下降，这样我们就可以一次性更新所有模型参数的猜测。

有多个参数要优化时，我们考虑损失表面，或者模型对一组可能的参数值的损失。

代码

import plotly.graph_objects as go

def mse_loss(theta, X, y_obs):
 y_hat = X @ theta
 return np.mean((y_hat - y_obs) ** 2) 

tips_with_bias = df.copy()
tips_with_bias["bias"] = 1
tips_with_bias = tips_with_bias[["bias", "total_bill"]]

uvalues = np.linspace(0, 2, 10)
vvalues = np.linspace(-0.1, 0.35, 10)
(u,v) = np.meshgrid(uvalues, vvalues)
thetas = np.vstack((u.flatten(),v.flatten()))

def mse_loss_single_arg(theta):
 return mse_loss(theta, tips_with_bias, df["tip"])

MSE = np.array([mse_loss_single_arg(t) for t in thetas.T])

loss_surface = go.Surface(x=u, y=v, z=np.reshape(MSE, u.shape))

ind = np.argmin(MSE)
optimal_point = go.Scatter3d(name = "Optimal Point",
 x = [thetas.T[ind,0]], y = [thetas.T[ind,1]], 
 z = [MSE[ind]],
 marker=dict(size=10, color="red"))

fig = go.Figure(data=[loss_surface, optimal_point])
fig.update_layout(scene = dict(
 xaxis_title = "theta0",
 yaxis_title = "theta1",
 zaxis_title = "MSE"), autosize=False, width=800, height=600)

fig.show()

我们还可以使用等高线图从上方可视化损失表面的鸟瞰图：

contour = go.Contour(x=u[0], y=v[:, 0], z=np.reshape(MSE, u.shape))
fig = go.Figure(contour)
fig.update_layout(
 xaxis_title = "theta0",
 yaxis_title = "theta1", autosize=False, width=800, height=600)

fig.show()

13.4.1 梯度向量

与以前一样，损失函数的导数告诉我们通往最小值的最佳方式。

在 2D（或更高维度）的表面上，向下（梯度）的最佳方式由向量描述。

数学旁注：偏导数 - 对于具有多个变量的方程，我们通过分别针对一个变量进行微分来进行偏导数。偏导数用$\partial$表示。直观地，我们想要看到函数在只改变一个变量的情况下如何变化，同时保持其他变量不变。- 以$f(x, y) = 3x^2 + y$为例，- 对 x 进行偏导数并将 y 视为常数，得到$\frac{\partial f}{\partial x} = 6x$ - 对 y 进行偏导数并将 x 视为常数，得到$\frac{\partial f}{\partial y} = 1$

对于参数值的向量$\vec{\theta} = \begin{bmatrix} \theta_{0} \\ \theta_{1} \\ \end{bmatrix}$，我们对每个参数的损失进行偏导数：$\frac{\partial L}{\partial \theta_0}$和$\frac{\partial L}{\partial \theta_1}$。

梯度向量因此为$\nabla_\theta L = \begin{bmatrix} \frac{\partial L}{\partial \theta_0} \\ \frac{\partial L}{\partial \theta_1} \\ \vdots \end{bmatrix}$其中$\nabla_\theta L$始终指向表面的下坡方向。

我们可以使用这个来更新我们的具有多个参数的 1D 梯度规则的模型。

• 回想一下我们的 1D 更新规则：$\theta^{(t+1)} = \theta^{(t)} - \alpha \frac{d}{d\theta}L(\theta^{(t)})$

• 对于具有多个参数的模型，我们使用向量来工作：$\begin{bmatrix} \theta_{0}^{(t+1)} \\ \theta_{1}^{(t+1)} \\ \vdots \end{bmatrix} = \begin{bmatrix} \theta_{0}^{(t)} \\ \theta_{1}^{(t)} \\ \vdots \end{bmatrix} - \alpha \begin{bmatrix} \frac{\partial L}{\partial \theta_{0}} \\ \frac{\partial L}{\partial \theta_{1}} \\ \vdots \\ \end{bmatrix}$

• 更紧凑地写成，$\vec{\theta}^{(t+1)} = \vec{\theta}^{(t)} - \alpha \nabla_{\vec{\theta}} L(\theta^{(t)})$

  • $\theta$ 是我们模型权重的向量

  • $L$ 是损失函数

  • $\alpha$ 是学习率（我们的是恒定的，但其他技术使用随时间减小的 $\alpha$）

  • $\vec{\theta}^{(t)}$ 是 $\theta$ 的当前值

  • $\vec{\theta}^{(t+1)}$ 是 $\theta$ 的下一个值

  • $\nabla_{\vec{\theta}} L(\theta^{(t)})$ 是在当前 $\vec{\theta}^{(t)}$ 处评估的损失函数的梯度

13.5 批量、小批量梯度下降和随机梯度下降

形式上，我们上面推导的算法称为批量梯度下降。对于算法的每次迭代，计算整个包含所有 $n$ 个数据点的批次的损失的导数。虽然这个更新规则在理论上效果很好，但在大多数情况下并不实用。对于大型数据集（可能有数十亿个数据点），在所有数据上找到梯度是非常耗费计算资源的；梯度下降会收敛缓慢，因为每次单独的更新都很慢。

小批量梯度下降试图解决这个问题。在小批量下降中，只使用子集数据来估计梯度。批量大小是每个子集中使用的数据点数。

每次完整“通过”数据称为训练周期。在小批量梯度下降的单个训练周期中，我们

• 计算第一个 x% 的数据的梯度。更新参数猜测。

• 计算下一个 x% 的数据的梯度。更新参数猜测。

• （省略）

• 计算最后 x% 的数据的梯度。更新参数猜测。

每个数据点在单个训练周期中只出现一次。然后我们进行多个训练周期，直到满意为止。

在最极端的情况下，我们可能选择只有 1 个数据点的批量大小——也就是说，每次更新步骤中只使用一个数据点来估计损失的梯度。这被称为随机梯度下降。在随机梯度下降的单个训练周期中，我们

• 计算第一个数据点的梯度。更新参数猜测。

• 计算下一个数据点的梯度。更新参数猜测。

• （省略）

• 计算最后一个数据点的梯度。更新参数猜测。

批量梯度下降是一种确定性技术——因为在每次更新迭代中都使用整个数据集，算法总是朝着损失曲面的最小值前进。相比之下，小批量和随机梯度下降都涉及一定的随机性。由于每次迭代更新参数 $\vec{\theta}$ 时只使用完整数据的子集，算法有可能不会在每次更新中朝着真正的损失最小值前进。从长远来看，这些随机技术仍然应该收敛到最优解。

下面的图表代表了从上方俯视的损失曲面。注意批量梯度下降直接朝向最优的 $\hat{\theta}$。相比之下，随机梯度下降在损失曲面上“跳跃”。这反映了每次更新步骤中抽样过程的随机性。

十四、Sklearn 和特征工程

原文：Sklearn and Feature Engineering

译者：飞龙

协议：CC BY-NC-SA 4.0

学习成果

• 应用sklearn库进行模型创建和训练

• 认识到特征工程作为提高模型性能的工具的价值

• 实现多项式特征生成和独热编码

• 了解模型复杂性、模型方差和训练误差之间的相互作用

到目前为止，我们已经对建模过程非常熟悉。我们介绍了损失的概念，用它来拟合多种类型的模型，并且最近扩展了我们的分析到多元回归。在这个过程中，我们通过推导出最佳模型参数的数学细节，艰难地走过了一段路。现在是时候让我们的生活变得更轻松一些了-让我们在代码中实现建模过程！

在本讲座中，我们将探讨两种模型拟合技术：

1. 将我们推导出的回归公式翻译成python

2. 使用python的sklearn包

有了我们手头的新编程框架，我们还将通过引入更复杂的特征来增强模型性能，为我们的模型增加复杂性。

14.1 在代码中实现推导出的公式

在本讲座中，我们将引用penguins数据集。

import pandas as pd
import seaborn as sns
import numpy as np

penguins = sns.load_dataset("penguins")
penguins = penguins[penguins["species"] == "Adelie"].dropna()
penguins.head()

	species	island	bill_length_mm	bill_depth_mm	flipper_length_mm	body_mass_g	sex
0	Adelie	Torgersen	39.1	18.7	181.0	3750.0	Male
1	Adelie	Torgersen	39.5	17.4	186.0	3800.0	Female
2	Adelie	Torgersen	40.3	18.0	195.0	3250.0	Female
4	Adelie	Torgersen	36.7	19.3	193.0	3450.0	Female
5	Adelie	Torgersen	39.3	20.6	190.0	3650.0	Male

我们的目标是预测给定企鹅的“flipper_length_mm”和“body_mass_g”时，“bill_depth_mm”的值。我们还将添加一个全为 1 的偏置列来表示我们模型的截距项。

# Add a bias column of all ones to `penguins`
penguins["bias"] = np.ones(len(penguins), dtype=int) 

# Define the design matrix, X...
X = penguins[["bias", "flipper_length_mm", "body_mass_g"]].to_numpy()

# ...as well as the target variable, y
Y = penguins[["bill_depth_mm"]].to_numpy()
# Converting X and Y to NumPy arrays avoids misinterpretation of column labels

在普通最小二乘法的讲座中，我们使用矩阵表示法表示多元线性回归。

$\hat{\mathbb{Y}} = \mathbb{X}\theta$

我们使用了几何方法来推导出最佳模型参数的以下表达式：

$\hat{\theta} = (\mathbb{X}^T \mathbb{X})^{-1}\mathbb{X}^T \mathbb{Y}$

这是大量的矩阵操作。我们如何在python中实现它呢？

这里有三个操作我们需要执行：矩阵相乘、求转置和求逆。

• 要进行矩阵乘法，使用@运算符

• 要进行转置，调用NumPy数组或DataFrame的.T属性

• 要计算逆矩阵，使用NumPy的内置方法np.linalg.inv

将这一切放在一起，我们可以计算存储在数组theta_hat中的最佳模型参数的 OLS 估计值。

theta_hat = np.linalg.inv(X.T @ X) @ X.T @ Y
theta_hat

array([[1.10029953e+01],
       [9.82848689e-03],
       [1.47749591e-03]])

要使用我们优化的参数值进行预测，我们将设计矩阵与参数向量进行矩阵乘法：

$\hat{\mathbb{Y}} = \mathbb{X}\theta$

Y_hat = X @ theta_hat
pd.DataFrame(Y_hat).head()

	0
0	18.322561
1	18.445578
2	17.721412
3	17.997254
4	18.263268

14.2 sklearn

通过将我们推导出的公式转化为代码，我们已经节省了大量时间（并避免了繁琐的计算）。但是，我们仍然需要自己编写线性代数的过程。

为了让生活变得更加轻松，我们可以转向sklearnpython库。sklearn是一个强大的机器学习工具库，在研究和工业中被广泛使用。它为我们提供了各种内置的建模框架和方法，因此在我们进行 Data 100 的过程中，我们将不断返回sklearn技术。

无论实现的具体模型类型是什么，sklearn都遵循一套标准的创建模型步骤。

1. 创建一个模型对象。这将生成模型类的一个新实例。你可以把它看作是对模型的标准“模板”进行新的“复制”。在代码中，这看起来像：

   my_model = ModelClass() 
   

2. 将模型拟合到X设计矩阵和Y目标向量。这将在不需要显式进行计算的情况下“在后台”计算出最佳的模型参数。然后，拟合的参数将存储在模型中以备将来进行预测使用：

   my_model.fit(X, Y) 
   

3. 使用拟合的模型对X输入数据进行预测，使用.predict。

   my_model.predict(X) 
   

要提取拟合的参数，我们可以使用：

my_model.coef_

my_model.intercept_ 

让我们在多元回归任务中将其付诸实践。

1. 初始化模型类的一个实例

sklearn存储了用于机器学习的有用模型的“模板”。我们通过制作一个这些模板的“副本”来开始建模过程以供我们自己使用。模型初始化看起来像ModelClass()，其中ModelClass是我们希望创建的模型类型。

现在，让我们使用LinearRegression()创建一个线性回归模型。

my_model现在是LinearRegression类的一个实例。你可以把它看作是线性回归模型的“想法”。我们还没有对它进行训练，所以它不知道任何模型参数，也不能用来进行预测。事实上，我们甚至还没有告诉它要用什么数据进行建模！它只是等待进一步的指示。

import sklearn.linear_model as lm

my_model = lm.LinearRegression()

2. 使用.fit训练模型

在模型可以进行预测之前，我们需要将其拟合到我们的训练数据中。当我们拟合模型时，sklearn将在后台运行梯度下降来确定最佳的模型参数。然后它会将这些模型参数保存到我们的模型实例中以备将来使用。

所有sklearn模型类都包括一个.fit方法，用于拟合模型。它接受两个输入：设计矩阵X和目标变量Y。

让我们从只有一个特征的模型开始：脚蹼长度。我们通过从DataFrame中提取"flipper_length_mm"列来创建一个设计矩阵X。

# .fit expects a 2D data design matrix, so we use double brackets to extract a DataFrame
X = penguins[["flipper_length_mm"]]
Y = penguins["bill_depth_mm"]

my_model.fit(X, Y)

LinearRegression()

**在 Jupyter 环境中，请重新运行此单元格以显示 HTML 表示或信任笔记本。

在 GitHub 上，HTML 表示无法呈现，请尝试使用 nbviewer.org 加载此页面。**

LinearRegression()

请注意，我们使用双括号来提取这一列。为什么使用双括号而不是单括号？.fit方法默认期望接收二维数据 - 一种包含行和列的数据。写penguins["flipper_length_mm"]会返回一个 1DSeries，导致sklearn出错。我们通过写penguins[["flipper_length_mm"]]来产生一个 2DDataFrame来避免这种情况。

我们的模型只用了三行代码就运行了梯度下降来确定最佳的模型参数！我们的单特征模型采用以下形式：

$\text{bill depth} = \theta_0 + \theta_1 \text{flipper length}$

请注意，LinearRegression将自动包括一个截距项。

拟合的模型参数被存储为模型实例的属性。my_model.intercept_将返回$\hat{\theta}_0$的值作为标量。my_model.coef_将以数组的形式返回所有值$\hat{\theta}_1, \hat{\theta}_1, ...$。因为我们的模型只包含一个特征，所以下面的单元格中只会看到$\hat{\theta}_1$的值。

# The intercept term, theta_0
my_model.intercept_

7.297305899612306

# All parameters theta_1, ..., theta_p
my_model.coef_

array([0.05812622])

3. 使用拟合的模型进行预测

现在模型已经训练好了，我们可以用它进行预测！为此，我们使用.predict方法。.predict接受一个参数：应该用来生成预测的设计矩阵。为了了解模型在训练集上的表现，我们会传入训练数据。或者，为了对未见过的数据进行预测，我们会传入一个未用于训练模型的新数据集。

在下面，我们调用.predict来在原始训练数据上生成模型预测。与之前一样，我们使用双括号来确保我们提取的是二维数据。

Y_hat_one_feature = my_model.predict(penguins[["flipper_length_mm"]])

print(f"The RMSE of the model is {np.sqrt(np.mean((Y-Y_hat_one_feature)**2))}")

The RMSE of the model is 1.1549363099239012

如果我们想要一个具有两个特征的模型呢？

$\text{bill depth} = \theta_0 + \theta_1 \text{flipper length} + \theta_2 \text{body mass}$

我们通过初始化一个新的模型对象，然后像以前一样调用.fit和.predict来重复这个三步过程。

# Step 1: initialize LinearRegression model
two_feature_model = lm.LinearRegression()

# Step 2: fit the model
X_two_features = penguins[["flipper_length_mm", "body_mass_g"]]
Y = penguins["bill_depth_mm"]

two_feature_model.fit(X_two_features, Y)

# Step 3: make predictions
Y_hat_two_features = two_feature_model.predict(X_two_features)

print(f"The RMSE of the model is {np.sqrt(np.mean((Y-Y_hat_two_features)**2))}")

The RMSE of the model is 0.9881331104079044

我们还可以看到，我们使用sklearn得到的预测与之前应用普通最小二乘公式时得到的预测相同！

代码

pd.DataFrame({"Y_hat from OLS":np.squeeze(Y_hat), "Y_hat from sklearn":Y_hat_two_features}).head()

	Y_hat from OLS	Y_hat from sklearn
0	18.322561	18.322561
1	18.445578	18.445578
2	17.721412	17.721412
3	17.997254	17.997254
4	18.263268	18.263268

14.3 特征工程

在课程的这一阶段，我们已经掌握了一些强大的技术来构建和优化模型。我们已经探讨了如何开发多变量模型，以及如何转换变量以帮助线性化数据集，并拟合这些模型以最大化它们的性能。

所有这些都是在一个主要的警告下完成的：到目前为止，我们所使用的回归模型都是输入变量的线性。我们假设我们的预测应该是一些线性变量的组合。虽然在某些情况下这很有效，但现实世界并不总是那么简单。我们将学习一种重要的方法来解决这个问题 - 特征工程 - 并考虑在这样做时可能出现的一些新问题。

特征工程是将原始特征转换为更具信息性的特征的过程，这些特征可以用于建模或 EDA 任务，并提高模型性能。

特征工程允许您：

• 捕获领域知识

• 使用线性模型表达非线性关系

• 在模型中使用非数值特征

14.4 特征函数

特征函数描述了我们对数据集中的原始特征应用的转换，以创建一个转换特征的设计矩阵。我们通常将特征函数表示为$\Phi$（想想：“phi”-ture 函数）。当我们将特征函数应用于我们的原始数据集$\mathbb{X}$时，结果$\Phi(\mathbb{X})$是一个经过转换的设计矩阵，可以用于建模。

例如，我们可以设计一个特征函数，计算现有特征的平方并将其添加到设计矩阵中。在这种情况下，我们现有的矩阵$[x]$被转换为$[x, x^2]$。其维度从 1 增加到 2。通常，特征化数据集的维度会增加，如此处所示。

特征函数引入的新特征然后可以用于建模。通常，我们使用符号$\phi_i$表示特征工程后的转换特征。

$\hat{y} = \theta_1 x + \theta_2 x^2$

$\hat{y}= \theta_1 \phi_1 + \theta_2 \phi_2$

在矩阵表示中，符号$\Phi$有时用于表示特征工程后的设计矩阵。请注意，在下面的用法中，$\Phi$现在是一个经过特征工程处理的矩阵，而不是一个函数。

$\hat{\mathbb{Y}} = \Phi \theta$

更正式地，我们将特征函数描述为将原始的$\mathbb{R}^{n \times p}$数据集$\mathbb{X}$转换为一个经过特征处理的$\mathbb{R}^{n \times p'}$数据集$\mathbb{\Phi}$的过程，其中$p'$通常大于$p$。

$\mathbb{X} \in \mathbb{R}^{n \times p} \longrightarrow \Phi \in \mathbb{R}^{n \times p'}$

14.5 独热编码

特征工程为设计更好的性能模型打开了一整套新的可能性。正如您将在实验室和家庭作业中看到的那样，特征工程是整个建模过程中最重要的部分之一。

特征工程的一个特别强大的用途是允许我们对非数值特征执行回归。独热编码是一种特征工程技术，它从分类数据生成数值特征，使我们能够使用通常的方法在数据上拟合回归模型。

为了说明这是如何工作的，我们将回顾以前讲座中的“小费”数据集。考虑数据集的“day”列：

代码

import numpy as np
tips = sns.load_dataset("tips")
tips.head()

	total_bill	tip	sex	smoker	day	time	size
0	16.99	1.01	Female	No	Sun	Dinner	2
1	10.34	1.66	Male	No	Sun	Dinner	3
2	21.01	3.50	Male	No	Sun	Dinner	3
3	23.68	3.31	Male	No	Sun	Dinner	2
4	24.59	3.61	Female	No	Sun	Dinner	4

乍一看，似乎不可能对这些数据拟合回归模型——我们无法直接对“太阳”进行任何数学运算。

为了解决这个问题，我们创建一个新表，其中包含原始“day”列中每个唯一值的特征。然后我们迭代“day”列。对于“day”中的每个条目，我们将新表中对应的特征填充为 1。所有其他特征都设置为 0。

sklearn的OneHotEncoder类（文档）提供了一种快速执行这种独热编码的方法。您将在实验室中详细探讨它的用法。现在，要认识到我们遵循的工作流程与我们使用LinearRegression类时非常相似：我们初始化一个OneHotEncoder对象，将其拟合到我们的数据，然后使用.transform来应用拟合的编码器。

from sklearn.preprocessing import OneHotEncoder

# Initialize a OneHotEncoder object
ohe = OneHotEncoder()

# Fit the encoder
ohe.fit(tips[["day"]])

# Use the encoder to transform the raw "day" feature
encoded_day = ohe.transform(tips[["day"]]).toarray()
encoded_day_df = pd.DataFrame(encoded_day, columns=ohe.get_feature_names_out())

encoded_day_df.head()

	day_Fri	day_Sat	day_Sun	day_Thur
0	0.0	0.0	1.0	0.0
1	0.0	0.0	1.0	0.0
2	0.0	0.0	1.0	0.0
3	0.0	0.0	1.0	0.0
4	0.0	0.0	1.0	0.0

然后，独热编码的特征可以用于设计矩阵来训练模型：

$\hat{y} = \theta_1 (\text{total}\textunderscore\text{bill}) + \theta_2 (\text{size}) + \theta_3 (\text{day}\textunderscore\text{Fri}) + \theta_4 (\text{day}\textunderscore\text{Sat}) + \theta_5 (\text{day}\textunderscore\text{Sun}) + \theta_6 (\text{day}\textunderscore\text{Thur})$

或者简写为：

$\hat{y} = \theta_1\phi_1 + \theta_2\phi_2 + \theta_3\phi_3 + \theta_4\phi_4 + \theta_5\phi_5 + \theta_6\phi_6$

现在，“day”特征（或者说，代表天的四个新布尔特征）可以用来拟合模型。

使用sklearn来拟合新模型，我们可以确定模型系数，从而了解每个特征对预测小费的影响。

from sklearn.linear_model import LinearRegression
data_w_ohe = tips[["total_bill", "size", "day"]].join(encoded_day_df).drop(columns = "day")
ohe_model = lm.LinearRegression(fit_intercept=False) #Tell sklearn to not add an additional bias column. Why?
ohe_model.fit(data_w_ohe, tips["tip"])

pd.DataFrame({"Feature":data_w_ohe.columns, "Model Coefficient":ohe_model.coef_})

	Feature	Model Coefficient
0	total_bill	0.092994
1	size	0.187132
2	day_Fri	0.745787
3	day_Sat	0.621129
4	day_Sun	0.732289
5	day_Thur	0.668294

例如，当查看day_Fri的系数时，我们可以了解星期五的事实对预测小费的影响有多大。

在独热编码时，要记住任何一组独热编码的列总是会加和为全为 1 的一列，表示偏置列。更正式地说，偏置列是 OHE 列的线性组合。

我们必须小心不要在我们的设计矩阵中包含这个偏置列。否则，模型中将存在线性依赖，这意味着$\mathbb{X}^T\mathbb{X}$将不再可逆，我们的 OLS 估计$\hat{\theta} = (\mathbb{X}^T\mathbb{X})^{-1}\mathbb{X}^T\mathbb{Y}$将失败。

为了解决这个问题，我们简单地省略了一个独热编码的列或不包括截距项。

任何一种方法都可以——我们仍然保留了与省略列相同的信息，即省略列是剩余列的线性组合。

14.6 多项式特征

我们已经遇到了几种情况，其中具有线性特征的模型在显示明显非线性曲率的数据集上表现不佳。

举个例子，考虑包含有关汽车信息的vehicles数据集。假设我们想要使用汽车的hp（马力）来预测其"mpg"（每加仑英里数的汽油里程）。如果我们可视化这两个变量之间的关系，我们会看到一个非线性的曲率。将线性模型拟合到这些变量会导致高（差）的 RMSE 值。

$\hat{y} = \theta_0 + \theta_1 (\text{hp})$

代码

pd.options.mode.chained_assignment = None 
vehicles = sns.load_dataset("mpg").dropna().rename(columns = {"horsepower": "hp"}).sort_values("hp")

X = vehicles[["hp"]]
Y = vehicles["mpg"]

hp_model = lm.LinearRegression()
hp_model.fit(X, Y)
hp_model_predictions = hp_model.predict(X)

import matplotlib.pyplot as plt

sns.scatterplot(data=vehicles, x="hp", y="mpg")
plt.plot(vehicles["hp"], hp_model_predictions, c="tab:red");

print(f"MSE of model with (hp) feature: {np.mean((Y-hp_model_predictions)**2)}")

MSE of model with (hp) feature: 23.943662938603108

为了捕捉数据集中的非线性，将非线性特征纳入其中是有意义的。让我们在回归模型中引入一个多项式项，$\text{hp}^2$。模型现在的形式是：

$\hat{y} = \theta_0 + \theta_1 (\text{hp}) + \theta_2 (\text{hp}^2)$

$\hat{y} = \theta_0 + \theta_1 \phi_1 + \theta_2 \phi_2$

我们如何拟合具有非线性特征的模型？我们可以使用与以前完全相同的技术：普通最小二乘法、梯度下降或sklearn。这是因为我们的新模型仍然是一个线性模型。尽管它包含非线性特征，但它在模型参数方面是线性的。我们以前所有拟合模型的工作都是在假设我们正在处理线性模型的情况下进行的。因为我们的新模型仍然是线性的，我们可以应用我们现有的方法来确定最佳参数。

# Add a hp^2 feature to the design matrix
X = vehicles[["hp"]]
X["hp^2"] = vehicles["hp"]**2

# Use sklearn to fit the model
hp2_model = lm.LinearRegression()
hp2_model.fit(X, Y)
hp2_model_predictions = hp2_model.predict(X)

sns.scatterplot(data=vehicles, x="hp", y="mpg")
plt.plot(vehicles["hp"], hp2_model_predictions, c="tab:red");

print(f"MSE of model with (hp^2) feature: {np.mean((Y-hp2_model_predictions)**2)}")

MSE of model with (hp^2) feature: 18.984768907617223

看起来好多了！通过引入一个平方特征，我们能够捕捉数据集的曲率。我们的模型现在是一个以我们的数据为中心的抛物线。请注意，相对于具有线性特征的原始模型，我们的新模型的误差已经减少了。

14.7 复杂性和过拟合

我们现在已经看到，特征工程使我们能够构建各种特征来提高模型的性能。特别是，我们看到设计更复杂的特征（在先前的vehicles数据中对hp进行平方）大大提高了模型捕捉非线性关系的能力。为了充分利用这一点，我们可能倾向于设计越来越复杂的特征。考虑以下三个不同阶数的模型（每个模型的最大指数幂）：

• 二次模型：$\hat{\text{mpg}} = \theta_0 + \theta_1 (\text{hp}) + \theta_2 (\text{hp}^2)$

• 三次模型：$\hat{\text{mpg}} = \theta_0 + \theta_1 (\text{hp}) + \theta_2 (\text{hp}^2) + \theta_3 (\text{hp}^3)$

• 四次模型：$\hat{\text{mpg}} = \theta_0 + \theta_1 (\text{hp}) + \theta_2 (\text{hp}^2) + \theta_3 (\text{hp}^3) + \theta_4 (\text{hp}^4)$

正如我们在上面的图中所看到的，随着每个额外的多项式项，均方误差继续减小。为了进一步可视化，让我们将模型从复杂度 0 增加到 6 进行绘制：

当我们使用我们的模型对用于拟合模型的相同数据进行预测时，我们发现随着每个额外的多项式项（随着我们的模型变得更复杂），MSE 会减少。训练误差是模型在生成来自用于训练目的的相同数据的预测时的误差。我们可以得出结论，随着模型复杂度的增加，训练误差会下降。

这看起来像是个好消息 - 在处理训练数据时，我们可以通过设计越来越复杂的模型来提高模型性能。

数学事实：给定$N$个重叠的数据点，我们总是可以找到一个通过所有这些点的$N-1$次多项式。

例如：总是存在一个 4 次多项式曲线，可以完美地模拟一个包含 5 个数据点的数据集！train_error

然而，高模型复杂性也带来了自己的一系列问题。在构建上述vehicles模型时，我们在整个数据集上训练了模型，然后评估了它们在同一数据集上的性能。实际上，我们很可能会在样本中训练模型，然后使用它对在训练期间未遇到的数据进行预测。

让我们通过一个更现实的例子来看看。假设我们有一个仅包含 6 个数据点的训练数据集，并希望训练一个模型，然后对不同的数据点进行预测。我们可能会倾向于制作一个非常复杂的模型（例如，5 次方），特别是考虑到它在左侧清晰地对训练数据进行了完美的预测。然而，如右侧图表所示，这个模型在整个数据集上的表现会非常糟糕！

上述现象被称为过拟合。该模型实际上只是记住了它在拟合时遇到的训练数据，导致它无法很好地对在训练期间未遇到的数据进行泛化。

此外，由于复杂模型对用于训练它们的特定数据集敏感，它们具有高方差。具有高方差的模型在训练不同数据集时往往会产生更大的变化。回到上面的例子，我们可以看到我们的 5 次方模型在拟合来自vehicles的不同 6 点样本时变化不稳定。

现在我们面临一个两难选择：我们知道我们可以通过增加模型复杂性来减少训练误差，但是过于复杂的模型开始过拟合，并且由于高方差无法重新应用于新的数据集。

我们可以看到，模型复杂性带来了明显的权衡。随着模型复杂性的增加，模型在训练数据上的误差减少。与此同时，模型的方差往往会增加。

这里的要点是：我们需要在模型的复杂性上取得平衡；我们希望模型能够泛化到“未见过”的数据。一个太简单的模型将无法捕捉我们感兴趣的变量之间的关键关系；一个太复杂的模型则有过拟合的风险。

这引出了一个问题：我们如何控制模型的复杂性？请关注我们的第 16 讲，交叉验证和正则化！

十五、人类背景和伦理的案例研究

原文：Case Study in Human Contexts and Ethics

译者：飞龙

协议：CC BY-NC-SA 4.0

学习成果

• 了解数据科学家面临的道德困境。

• 了解如何使用有关数据的上下文知识批判模型。

免责声明：以下章节讨论了结构性种族主义的问题。本章中的一些内容可能比较敏感，可能或可能不是收集材料的学生的意见、想法和信念。Data 100 课程工作人员尽最大努力只呈现与教授课程相关的信息。

**注意：**由于讲座中提出的一些论点的微妙性质，强烈建议您观看讲座录像，以便充分参与和理解材料。课程笔记将具有相同的更广泛结构，但绝不是全面的。

让我们沉浸在一个名为库克县评估员办公室（CCAO）的组织的数据科学家的现实故事中。他们的工作是估算房屋的价值以分配财产税。这是因为该地区的税收负担是由房屋的估算价值决定的，这与其价格不同。由于价值随时间变化，且没有明显的价值指标，他们创建了一个模型来估算房屋的价值。在本章中，我们将深入探讨偏见模型存在的问题，对人类生活的后果，以及我们如何从这个例子中学习以做得更好。

15.1 问题

《芝加哥论坛报》的一份报告揭露了一个重大丑闻：该团队表明该模型延续了一个高度累进的税收体系，不成比例地加重了库克县的非裔美国人和拉丁裔房主的负担。他们是如何知道的呢？

在住房评估领域，评估员使用的标准指标是：离散系数和价格相关差异。这些指标已经在该领域的专家进行了严格测试，超出了我们课程的范围。对库克县价格计算这些指标显示，CCAO 制定的价格不在可接受范围内（见上图）。这本身并不是故事的结束，但是一个很好的指示有些不对劲。

这促使他们调查模型本身是否产生了公平的税率。显然，当考虑到房主的收入时，他们发现该模型实际上产生了累进的税率（见上图）。如果百分比税率对收入较低的个人更高，则税率是累进的。如果百分比税率对收入较高的个人更高，则税率是累进的。

进一步的调查表明，这个系统不仅对收入的轴线上的人不公平，也对种族的轴线上的人不公平（见上图）。一栋房产被低估或高估的可能性高度依赖于业主的种族，这让许多房主感到不安。

15.1.1 重点：上诉

这到底是什么导致了这种情况？一个全面的答案超出了模型。归根结底，这些都是有很多活动部分的真实系统。其中之一是上诉系统。房主会收到 CCAO 评估的房屋价值的邮件，房主可以选择向一组选举官员上诉，试图改变他们的房屋价值清单，从而改变他们被征税的金额。从理论上讲，这听起来是一个非常公平的系统：有人监督房屋的最终定价，而不仅仅是一个算法。然而，它最终加剧了问题。

“上诉是一件好事，”估价和上诉副主管托马斯·雅科内蒂在一次采访中说。“这里的目标是公平。我们制定了这些数字。我们可以改变它们。”

在这里，我们可以从批判种族理论中汲取教训。表面上，每个人都有法定权利尝试上诉是不可否认的。然而，并非每个人都有平等的能力这样做。那些有钱雇佣税务律师为他们上诉的人尝试和成功的机会大大提高（见上图）。这种模式是潜在腐败的深层制度模式的一部分。

上诉的房主相对于没有上诉的房主通常被低估（见上图）。收入较高的人支付较少的房产税，税务律师能够因为他们在上诉中的角色而发展业务，政客通常与上述税务律师和富裕的房主有社会联系。所有这些利益相关者都有理由宣传这种模式作为公平制度的一个组成部分。在这里提出问题是有价值的：一个表面上看起来公平的系统在仔细观察后实际上可能是不公平的。

15.1.2 人类影响

住房模式的影响超出了房屋所有权和税收的范围。歧视性做法在美国有着悠久的历史，而这种模式则服务于延续这一事实。直到今天，芝加哥是美国最种族隔离的城市之一（来源）。这些因素对我们作为数据科学家来说至关重要。

15.1.3 聚焦：房地产和种族的交汇

住房一直是美国历史上种族不平等的持续来源，也是其他因素之一。它是不平等产生和再生产的主要领域之一。起初，吉姆·克劳法明确禁止有色人种进入学校、公共设施等。

今天，尽管民权方面取得了进展，但美国许多地方的法律精神仍然存在。房地产行业在 20 世纪 20 年代和 30 年代被“专业化”，渴望成为一门由严格方法和原则指导的科学，如下所述：

• 红线政策：使在特定被编码为“高风险”（红色）的社区购买联邦支持抵押贷款变得困难或不可能。

  • 根据这些制造商的说法，是什么使它们“高风险”。

  • 种族隔离不仅是联邦政策的结果，而且是房地产专业人士发展起来的。

• 方法集中在创建客观的评级系统（信息技术）来评估房产价值，这些系统将种族编码为估值因素（见下图），

  • 这反过来影响了联邦政策和实践。

来源：Colin Koopman，《我们如何成为我们的数据》（2019 年）第 137 页

15.2 回应：库克县开放数据倡议

回应始于政治。新的评估员弗里茨·凯吉当选并制定了两个目标的新任务：

1. 财产税的分配公平，意味着同等价值的财产在评估过程中受到同等对待。

2. 创建一个新的数据科学办公室。

15.2.1 问题/问题的形成

驱动问题

• 我们想知道什么？

• 我们试图解决什么问题？

• 我们想要测试什么假设？

• 我们的成功指标是什么？

数据科学办公室通过重新定义他们的目标来开始。

1. 准确、一致和公正地评估房屋价值

   • 遵循国际标准（离散系数）

   • 尽可能准确地预测所有房屋的价值

2. 创建一个强大的管道，能够准确评估规模化的财产价值，并且公平

   • 打破腐败的循环（评审委员会的上诉程序）

   • 消除累退性

   • 在所有利益相关者中建立对系统的信任

定义：公平和透明度

库克县评估员办公室给出的定义如下：

• 公平：我们的管道能够准确评估财产价值，考虑到地理、信息等方面的差异。

• 透明度：数据科学部门分享和解释管道结果和决策的能力，向内部和外部利益相关者

请注意办公室如何以准确性来定义“公平”。因此，问题——使系统更公平——已经以数据科学家易于处理的术语来表述：使评估更准确。

这里的想法是，如果模型更准确，它也会（或许必然会）变得更公平，这是一个很大的假设。从某种意义上说，存在两个不同的问题——做出准确的评估，建立一个公平的系统。

目标的定义方式导致我们提出了一个问题：准确评估财产价值实际上意味着什么，以及“规模”在其中扮演了什么角色？

1. 房屋价值的评估是什么？

2. 一项评估比另一项更准确的原因是什么？

3. 一批评估比另一批更准确的原因是什么？

以上每个问题都引发了一系列更多的问题。仅考虑第一个问题，一个答案可能是评估是对房屋价值的估计。这带来了更多的疑问：房屋的价值是多少？是什么决定了它？我们怎么知道？在这个课程中，我们认为它是房屋的市场价值。

15.2.2 数据获取和清理

驱动问题

• 我们有什么数据，我们需要什么数据？

• 我们将如何取样更多的数据？

• 我们的数据是否代表我们想研究的人口？

数据科学家还对他们最初的销售数据进行了批判性审查：

并提出了以下问题：

1. 这些数据是如何收集的？

2. 这些数据是何时收集的？

3. 谁收集了这些数据？

4. 数据收集的目的是什么？

5. 特定类别是如何创建的，为什么创建的？

例如，属性在数据中出现的可能性不同，芝加哥洪水平原地理区域的住房数据比其他地区少。

这些特征甚至可以以不同的速度报告。改善房屋，往往会增加房产价值，但业主们不太可能报告这一点。

此外，他们发现低收入社区的缺失数据更多。

15.2.3 探索性数据分析

驱动问题

• 我们的数据是如何组织的，它包含了什么？

• 我们已经有相关的数据了吗？

• 数据中存在什么偏见、异常或其他问题？

• 我们如何转换数据以实现有效的分析？

在建模步骤之前，他们调查了许多关键问题：

1. 哪些属性对销售价格的预测最有帮助？

2. 数据是否均匀分布？

3. 所有社区的数据都是最新的吗？所有社区的粒度都一样吗？

4. 一些社区的数据是否缺失或过时？

首先，他们发现某些特征的影响，比如卧室数量，在确定某些社区内房屋价值方面比其他社区更有影响力。这让他们知道应该根据社区使用不同的模型。

他们还注意到低收入社区的数据不足。这让他们知道他们需要开发新的数据收集实践，包括寻找新的数据来源。

15.2.4 预测和推断

驱动问题

• 数据对世界有何说法？

• 它是否回答了我们的问题或准确解决了问题？

• 我们的结论有多可靠，我们能相信这些预测吗？

CCA0 不是使用单一模型来预测未售出房产的销售价格（“公平市场价”），而是使用机器学习模型来发现使用已知销售价格和相似和附近房产特征的模式。它为每个乡镇使用不同的模型权重。

与传统的大规模评估相比，CCA0 的新方法更加细致，并对社区变化更加敏感。

在这里，我们可能会问为什么任何特定的个人应该相信该模型对他们的房产是准确的？

这让我们认识到，CCA0 依赖于其“透明度”的表现（将数据、模型、管道放到 GitLab 上）来促进公众的信任，这将有助于将“准确评估”的产出与“公平”相提并论。

在评估我们的模型时，准确性、公平性和我们倾向使用的指标之间的关系还有很多需要讨论的地方。鉴于论点的微妙性质，建议您查看相应的讲座，因为课程笔记对于讲座的这一部分并不那么全面。

15.2.5 报告决策和结论

驱动问题

• 系统对每个目标有多成功？

  • 模型的准确性/一致性

  • 消除回归性和建立信任的公平和透明度

• 你怎么知道的？

模型并不是终点。新办公室仍然向房主发送他们的房屋评估报告，但现在他们会考虑从房主那里得到的数据。办公室本身正在撰写更详细的报告，以使信息民主化。市政厅和其他面向公众的宣传活动有助于让整个社区参与到住房评估的过程中，而不是仅限于少数人参与。****

15.3 主要收获

1. 准确性是公平系统的必要条件，但不是充分条件。

2. 公平和透明度是依赖于环境和社会技术概念的。

3. 学会与环境一起工作，并考虑你的数据分析将如何重塑它们。

4. 牢记数据分析的力量和局限性。

15.4 数据科学实践的教训

1. 问题/问题的制定

   • 谁负责框定问题？

   • 谁是利益相关者？他们如何参与问题的框架？

   • 你能提供什么？你的立场如何影响你对问题的理解？

   • 你正在利用哪些叙事？

2. 数据获取和清理

   • 数据从哪里来？

   • 谁收集了它？为了什么目的？

   • 使用了什么样的收集和记录系统和技术？

   • 过去如何使用这些数据？

   • 对于数据访问有什么限制，以及是什么使您能够访问？

3. 探索性数据分析和可视化

   • 在这些数据中，有哪些个人或群体身份变得突出？

   • 哪些变量变得突出，它们之间有什么样的关系？

   • 任何可见的关系是否会导致对特定社区可能有害的争论？

4. 预测和推断

   • 预测或推断在世界上起到什么作用？

   • 结果对预期目的有用吗？

   • 是否有基准可以比较结果？

   • 你的预测和推断如何依赖于模型所在的更大系统？

5. 报告、决策和解决方案

   • 我们如何知道我们是否实现了我们的目标？

   • 你的工作如何融入更广泛的文献？

   • 你的工作在哪些方面与现状一致或不一致？

   • 你的结论是否合理？