《机器学习的数学基础》第9章"线性回归"

第 9 章 线性回归

Linear Regression

在本章中，我们将应用第 2、5、6 和 7 章中的数学概念来解决线性回归（曲线拟合 curve fitting）问题。在回归中，我们的目标是找到一个函数 \(f\)，它将输入 \(x \in \mathbb{R}^D\) 映射到相应的函数值 \(f(x) \in \mathbb{R}\)。

我们假设给定了一组训练输入 \(x_n\) 及其对应的带噪观测值 \(y_n = f(x_n) + \varepsilon\)，其中 \(\varepsilon\) 是独立同分布（i.i.d.）的随机变量，用来描述测量/观测噪声以及可能未建模的过程（本章中我们不再考虑后者）。在整个章节中，我们假设噪声是零均值高斯噪声。

我们的任务是找到一个函数，该函数不仅能够拟合训练数据，还能够很好地推广到预测训练数据之外输入位置的函数值（见第 8 章）。图 9.1 展示了这样一个回归问题的例子。典型的回归场景如图 9.1(a)：对于一些输入值 \(x_n\)，我们观测到带噪声的函数值 \(y_n = f(x_n) + \varepsilon\)。任务是推断出生成这些数据的函数 \(f\)，并且能很好地推广到新的输入位置。图 9.1(b) 给出了一个可能的解答，其中我们同时展示了以函数值 \(f(x)\) 为中心的三个分布，用来表示数据中的噪声。

回归是机器学习中的一个基本问题，回归问题广泛出现在各种研究领域和应用中，包括：

• 时间序列分析（如系统辨识）
• 控制与机器人（如强化学习、前向/逆向模型学习）
• 优化（如线搜索、全局优化）
• 深度学习应用（如电子游戏、语音转文本转换、图像识别、自动视频标注）

此外，回归还是分类算法中的关键组成部分。寻找回归函数需要解决多种问题，包括：

1. 模型选择与参数化
   • 给定一个数据集，哪些函数类别（如多项式）是建模的良好候选？
   • 具体的参数化（如多项式的次数）应该如何选择？
   • 第 8.6 节讨论的模型选择方法，可以帮助我们比较不同模型，从中找到能够合理解释训练数据的最简单模型。
2. 寻找合适的参数
   • 在选定回归函数模型后，如何找到合适的模型参数？
   • 我们需要研究不同的损失/目标函数（它们决定了“拟合得好”的含义），以及能够最小化该损失的优化算法。
3. 过拟合与模型选择
   • 当回归函数对训练数据“拟合得过好”但无法推广到未见过的测试数据时，就会出现过拟合。
   • 过拟合通常出现在底层模型（或其参数化）过于灵活、表达能力过强的情况下（见第 8.6 节）。
   • 我们将探讨过拟合的根本原因，并在线性回归的背景下讨论减轻过拟合的方法。
4. 损失函数与参数先验的关系
   • 损失函数（优化目标）通常由概率模型引导并激发。
   • 我们将探讨损失函数与其背后引入这些损失的先验假设之间的联系。
5. 不确定性建模
   • 在实际问题中，我们只能获得有限（尽管可能较大）数量的训练数据，用于选择模型类别和相应参数。
   • 由于这些有限的训练数据无法覆盖所有可能情况，我们可能需要刻画剩余的参数不确定性，以便在测试时得到模型预测的置信度度量。
   • 训练集越小，不确定性建模就越重要。
   • 一致的不确定性建模能够为模型预测提供置信区间。

在接下来的内容中，我们将使用第 3、5、6 和 7 章中的数学工具来解决线性回归问题。我们将讨论极大似然估计（MLE）和最大后验估计（MAP）以找到最优模型参数。利用这些参数估计，我们将简要考察泛化误差和过拟合。在本章的最后，我们将讨论贝叶斯线性回归，它能在更高层次上推理模型参数，从而避免极大似然和最大后验估计中遇到的一些问题。

9.1 问题表述

Problem Formulation

由于存在观测噪声，我们将采用一种概率方法，并通过似然函数显式地对噪声进行建模。更具体地说，在本章中，我们考虑一个回归问题，其似然函数为: \[ p(y \mid x) = \mathcal{N}\big(y \mid f(x), \sigma^2 \big). \tag{9.1} \] 这里，\(x \in \mathbb{R}^D\) 是输入，\(y \in \mathbb{R}\) 是带有噪声的函数值（目标）。根据式 (9.1)，\(x\) 与 \(y\) 之间的函数关系为

\[ y = f(x) + \varepsilon , \tag{9.2} \] 其中 \(\varepsilon \sim \mathcal{N}(0, \sigma^2)\) 是独立同分布 (i.i.d.) 的高斯测量噪声，均值为 0，方差为 \(\sigma^2\)。

个人注：

严格来说，似然函数（likelihood function） 是指在给定参数 \(\theta\) 的情况下，观测数据出现的概率。因此形式应该写成 \[ p(y \mid x, \theta) = \mathcal{N}\big(y \mid f(x;\theta), \sigma^2 \big), \tag{9.1} \] 其中 \(\theta\) 表示模型参数（比如线性回归中的权重向量和偏置，或者神经网络中的权重）。在很多教材（尤其是深度学习、贝叶斯建模的入门部分），作者会 省略参数 \(\theta\)，默认 \(f(x)\) 依赖于 \(\theta\)。这是出于简化符号的考虑，因为 \(\theta\) 总是“隐含存在”的。

我们的目标是找到一个函数，它既接近（类似于）生成数据的未知函数 \(f\)，又具有良好的泛化能力。

在本章中，我们专注于参数化模型。也就是说，我们选择一个参数化的函数，并寻找合适的参数 \(\theta\)，使其能够很好地建模数据。暂时我们假设噪声方差 \(\sigma^2\) 已知，并专注于学习模型参数 \(\theta\)。在线性回归中，我们考虑参数 \(\theta\) 以线性形式出现在模型中的特殊情况。线性回归的一个例子为：

\[ p(y \mid x, \theta) = \mathcal{N}\big(y \mid x^{\top}\theta, \sigma^2 \big), \tag{9.3} \]

$$ \[\begin{equation} \;\;\Longleftrightarrow\;\; y = x^{\top}\theta + \varepsilon, \quad \varepsilon \sim \mathcal{N}(0, \sigma^{2}). \tag{9.4} \end{equation}\] $$

其中 \(\theta \in \mathbb{R}^D\) 是我们需要估计的参数。由式 (9.4) 描述的函数类是穿过原点的直线。在式 (9.4) 中，我们选择了如下参数化形式： \[ f(x) = x^{\top}\theta. \] 式 (9.3) 中的似然函数是 \(y\) 的概率密度函数，在 \(x^{\top}\theta\) 处进行评估。需要注意的是，不确定性的唯一来源是观测噪声（因为在式 (9.3) 中，\(x\) 和 \(\theta\) 被假设为已知）。如果没有观测噪声，\(x\) 与 \(y\) 之间的关系将是确定性的，而式 (9.3) 将退化为一个狄拉克 delta 函数。

F9.2

例 9.1 对于 \(x, \theta \in \mathbb{R}\)，线性回归模型 (9.4) 描述的是直线（线性函数），参数 \(\theta\) 就是直线的斜率。图 9.2(a) 展示了不同参数取值下的一些示例函数。

线性回归模型 (9.3)–(9.4) 不仅在参数上是线性的，在输入 \(x\) 上也是线性的。图 9.2(a) 展示了此类函数的例子。我们将在后面看到，对于非线性变换 \(\phi\)，模型 \[ y = \phi^{\top}(x)\theta \] 依然是一个线性回归模型，因为“线性回归”指的是在参数上是线性的模型，即模型通过输入特征的线性组合来描述函数。这里的“特征”是输入 \(x\) 的一种表示 \(\phi(x)\)。在接下来的内容中，我们将更详细地讨论如何找到合适的参数 \(\theta\)，以及如何评估某组参数是否“效果良好”。在目前的讨论中，我们假设噪声方差 \(\sigma^2\) 是已知的。

9.2 参数估计

考虑线性回归的设定 (9.4)，假设我们给定一个训练集\(D := \{(x_1, y_1), \ldots, (x_N, y_N)\}\)其中包含 \(N\) 个输入 \(x_n \in \mathbb{R}^D\) 以及对应的观测/目标 \(y_n \in \mathbb{R}, \; n = 1, \ldots, N\)。对应的图模型如图 9.3 所示。注意到，给定各自的输入 \(x_i, x_j\) 时，\(y_i\) 和 \(y_j\) 条件独立，因此似然可以因子化为: \[ \begin{align} p(Y \mid X, \theta) & = p(y_{1}, \ldots, y_{N} \mid x_{1}, \ldots, x_{N}, \theta) \tag{9.5a} \\ & = \prod_{n=1}^{N} p(y_{n} \mid x_{n}, \theta) \\ & = \prod_{n=1}^{N} \mathcal{N}\!\left(y_{n} \mid x_{n}^{\top}\theta, \sigma^{2}\right) \tag{9.5b} \end{align} \] 其中我们定义 \(X := \{x_1, \ldots, x_N\}\)，\(Y := \{y_1, \ldots, y_N\}\) 分别为训练输入集合和对应的目标集合。由于噪声分布的存在，似然以及各个因子 \(p(y_n \mid x_n, \theta)\) 都是高斯分布；见 (9.3)。

接下来，我们将讨论如何为线性回归模型 (9.4) 找到最优参数 \(\theta^* \in \mathbb{R}^D\)。一旦找到参数 \(\theta^*\)，我们就可以利用该参数估计来进行预测：对于任意测试输入 \(x^*\)，其对应目标 \(y^*\) 的分布为 \[ p(y^* \mid x^*, \theta^*) = \mathcal{N}\!\big(y^* \mid (x^*)^\top \theta^*, \sigma^2\big). \tag{9.6} \] 在下文中，我们将考察通过最大化似然来进行参数估计的方法，这是我们在第 8.3 节中已经部分涉及过的一个主题。

9.2.1 极大似然估计

Maximum Likelihood Estimation

一种广泛使用的求解期望参数 \(\theta_{\text{ML}}\) 的方法是 极大似然估计，其思想是找到能最大化似然函数 (9.5b) 的参数 \(\theta_{\text{ML}}\)。直观上，最大化似然就意味着：在给定模型参数的条件下，使训练数据的预测分布最大。最大似然参数可以写为： \[ \theta_{\text{ML}} = \arg\max_{\theta} p(Y \mid X, \theta). \tag{9.7} \]

备注：似然函数 \(p(y \mid x, \theta)\) 并不是关于 \(\theta\) 的概率分布：它只是参数 \(\theta\) 的一个函数，但并不一定积分为 1（即它是未归一化的），甚至可能不可积。然而，在 (9.7) 中的似然函数，是关于 \(y\) 的一个归一化的概率分布。

为了找到能最大化似然的参数 \(\theta_{\text{ML}}\)，我们通常会进行 梯度上升（或者对负似然做 梯度下降）。但是，在这里讨论的线性回归情形下，存在一个 闭式解，因此无需迭代的梯度下降。实际应用中，我们通常不会直接最大化似然，而是对似然函数取对数，并最小化 负对数似然。

备注（对数变换）：由于似然函数 (9.5b) 是 \(N\) 个高斯分布的乘积，因此取对数变换非常有用，因为： (a) 它避免了数值下溢的问题； (b) 它使得求导规则更简单。

更具体地说，当我们相乘 \(N\) 个概率时（\(N\) 是数据点的个数），会出现数值下溢的问题，因为我们无法表示极小的数，例如 \(10^{-256}\)。此外，对数变换能把乘积转化为对数概率的和，这样相应的梯度就是单个梯度的和，而不必反复应用乘积法则 (5.46) 去计算 \(N\) 项乘积的梯度。

为了找到我们线性回归问题的最优参数 \(\theta_{\text{ML}}\)，我们最小化 负对数似然： \[ - \log p(Y \mid X, \theta) = - \log \prod_{n=1}^{N} p(y_n \mid x_n, \theta) = - \sum_{n=1}^{N} \log p(y_n \mid x_n, \theta), \tag{9.8} \] 其中我们利用了在训练集上独立性假设，使得似然函数 (9.5b) 可以在数据点数目上进行因子分解。在 线性回归模型 (9.4) 中，由于高斯加性噪声项的存在，似然函数是高斯分布，因此我们得到：

\[ \log p(y_n \mid x_n, \theta) = -\frac{1}{2\sigma^2} (y_n - x_n^\top \theta)^2 + \text{const}, \tag{9.9} \] 其中常数项包括了所有与 \(\theta\) 无关的部分。将 (9.9) 代入负对数似然 (9.8)，并忽略常数项后，我们得到：

\[ \begin{align} L(\theta) &:= \frac{1}{2\sigma^{2}} \sum_{n=1}^{N} \left( y_{n} - x_{n}^{\top}\theta \right)^{2} \tag{9.10a} \\ &= \frac{1}{2\sigma^{2}} (y - X\theta)^{\top}(y - X\theta) = \frac{1}{2\sigma^{2}} \, \lVert y - X\theta \rVert^{2}. \tag{9.10b} \end{align} \]

其中我们定义 设计矩阵（design matrix） \[ X := [x_1, \dots, x_N]^\top \in \mathbb{R}^{N \times D} \] 为训练输入的集合，且 \[ y := [y_1, \dots, y_N]^\top \in \mathbb{R}^N \] 为收集所有训练目标的向量。注意，设计矩阵 \(X\) 的第 \(n\) 行对应于训练输入 \(x_n\)。在公式 (9.10b) 中，我们利用了这样一个事实：观测值 \(y_n\) 与对应模型预测 \(x_n^\top \theta\) 之间的平方误差之和，等于向量 \(y\) 与 \(X\theta\) 之间的平方距离。

根据公式 (9.10b)，我们现在得到了需要优化的负对数似然函数的具体形式。我们立刻可以看到，(9.10b) 是关于参数 \(\theta\) 的二次函数。这意味着我们可以找到唯一的全局解 \(\theta_{\text{ML}}\)，从而最小化负对数似然 \(L\)。找到全局最优解的方法是：计算 \(L\) 对 \(\theta\) 的梯度，将其设为 0，然后解出 \(\theta\)。

利用第 5 章的结果，我们计算 \(L\) 关于参数的梯度如下： \[ \begin{align} \frac{dL}{d\theta} &= \frac{d}{d\theta} \left[ \frac{1}{2\sigma^{2}} (y - X\theta)^{\top}(y - X\theta) \right] \tag{9.11a} \\[6pt] &= \frac{1}{2\sigma^{2}} \frac{d}{d\theta} \left( y^{\top}y - 2y^{\top}X\theta + \theta^{\top}X^{\top}X\theta \right) \tag{9.11b} \\[6pt] &= \frac{1}{\sigma^{2}} \left( -y^{\top}X + \theta^{\top}X^{\top}X \right) \in \mathbb{R}^{1 \times D}. \tag{9.11c} \end{align} \] 极大似然估计量 \(\theta_{\text{ML}}\) 通过解方程 \(\frac{dL}{d\theta} = 0^\top\)（必要的最优条件）得到： \[ \begin{align} \frac{dL}{d\theta} = 0^{\top} &\;\;\;\;\;\;\;\;\; \Longleftrightarrow \;\;\;\;\;\;\;\;\; \theta^{\top}_{\mathrm{ML}} X^{\top}X = y^{\top}X \tag{9.12a} \\[6pt] &\;\;\;\;\;\;\;\;\; \Longleftrightarrow \;\;\;\;\;\;\;\;\; \theta^{\top}_{\mathrm{ML}} = y^{\top}X (X^{\top}X)^{-1} \tag{9.12b} \\[6pt] &\;\;\;\;\;\;\;\;\; \Longleftrightarrow \;\;\;\;\;\;\;\;\; \theta_{\mathrm{ML}} = (X^{\top}X)^{-1} X^{\top} y . \tag{9.12c} \end{align} \]

我们之所以能在 (9.12a) 的等式右乘 \((X^\top X)^{-1}\)，是因为当矩阵 \(X\) 的秩 \(\mathrm{rk}(X) = D\) 时，\(X^\top X\) 是正定的。

备注 1：将梯度设为 \(0^\top\) 是充分必要条件，并且我们得到的是全局最小值，因为 Hessian 矩阵 \[ \nabla^2_\theta L(\theta) = X^\top X \in \mathbb{R}^{D \times D} \] 是正定的。

备注 2：最大似然解 (9.12c) 要求我们解一个形如 \(A\theta = b\) 的线性方程组，其中 \[ A = X^\top X, \quad b = X^\top y. \]

例 9.2（拟合直线） 我们来看图 9.2，其中我们的目标是用极大似然估计将一条直线\(f(x) = \theta x\)（\(\theta\) 是未知的斜率）拟合到一个数据集上。该模型类别中的一些函数（直线）展示在图 9.2(a) 中。对于图 9.2(b) 所示的数据集，我们利用公式 (9.12c) 得到斜率参数 \(\theta\) 的极大似然估计，并在图 9.2(c) 中得到对应的最大似然直线函数。

带特征的极大似然估计

Maximum Likelihood Estimation with Features

到目前为止，我们考虑的是 (9.4) 中描述的线性回归设定，它允许我们用极大似然估计将直线拟合到数据上。然而，在拟合更复杂的数据时，直线的表达能力是不够的。幸运的是，线性回归在其框架下为我们提供了一种方法，可以拟合非线性函数：

由于“线性回归”仅仅指的是“对参数是线性的”，我们可以对输入 \(x\) 进行任意的非线性变换 \(\phi(x)\)，然后再对这个变换后的分量进行线性组合。相应的线性回归模型为 \[ p(y \mid x, \theta) = \mathcal{N}\big(y \mid \phi^\top(x)\theta, \sigma^2 \big) \quad \Longleftrightarrow \quad y = \phi^\top(x)\theta + \varepsilon = \sum_{k=0}^{K-1} \theta_k \phi_k(x) + \varepsilon , \tag{9.13} \] 其中 \(\phi : \mathbb{R}^D \to \mathbb{R}^K\) 是输入 \(x\) 的一个（非线性）变换，而 \(\phi_k : \mathbb{R}^D \to \mathbb{R}\) 是特征向量 \(\phi\) 的第 \(k\) 个分量。需要注意的是，模型参数 \(\theta\) 仍然只以线性形式出现。

例 9.3（多项式回归） 我们关心的是一个回归问题： \[ y = \phi^\top(x)\theta + \varepsilon, \] 其中 \(x \in \mathbb{R}\)，\(\theta \in \mathbb{R}^K\)。在这种情况下常用的一种变换是 \[ \phi(x) = \begin{bmatrix} \phi_0(x) \\ \phi_1(x) \\ \vdots \\ \phi_{K-1}(x) \end{bmatrix} = \begin{bmatrix} 1 \\ x \\ x^2 \\ x^3 \\ \vdots \\ x^{K-1} \end{bmatrix} \in \mathbb{R}^K. \tag{9.14} \] 这意味着我们将原本一维的输入空间“提升”到一个 \(K\) 维特征空间，该空间由所有单项式 \(x^k\) （其中 \(k = 0, \ldots, K-1\)）组成。利用这些特征，我们可以在线性回归框架下建模次数不超过 \(K-1\) 的多项式：

一个 \(K-1\) 次多项式是 \[ f(x) = \sum_{k=0}^{K-1} \theta_k x^k = \phi^\top(x)\theta, \tag{9.15} \] 其中 \(\phi\) 定义如 (9.14)，\(\theta = [\theta_0, \ldots, \theta_{K-1}]^\top \in \mathbb{R}^K\) 包含了参数 \(\theta_k\)（线性参数）。

现在我们来看看线性回归模型 (9.13) 中参数 \(\theta\) 的最大似然估计。我们考虑训练输入 \(x_n \in \mathbb{R}^D\) 和目标 \(y_n \in \mathbb{R}, \; n = 1, \dots, N\)，并定义特征矩阵（设计矩阵）为 \[ \Phi := \begin{bmatrix} \phi^{\top}(x_1) \\ \vdots \\ \phi^{\top}(x_N) \end{bmatrix} = \begin{bmatrix} \phi_0(x_1) & \cdots & \phi_{K-1}(x_1) \\ \phi_0(x_2) & \cdots & \phi_{K-1}(x_2) \\ \vdots & \ddots & \vdots \\ \phi_0(x_N) & \cdots & \phi_{K-1}(x_N) \end{bmatrix} \in \mathbb{R}^{N \times K}, \tag{9.16} \] 其中 \(\Phi_{ij} = \phi_j(x_i)\)，并且 \(\phi_j : \mathbb{R}^D \to \mathbb{R}\)。

例 9.4 （二次多项式的特征矩阵） 对于二次多项式以及训练点 \(x_n \in \mathbb{R}, \; n=1,\dots,N\)，其特征矩阵为 \[ \Phi = \begin{bmatrix} 1 & x_1 & x_1^2 \\ 1 & x_2 & x_2^2 \\ \vdots & \vdots & \vdots \\ 1 & x_N & x_N^2 \end{bmatrix}. \tag{9.17} \]

在 (9.16) 定义的特征矩阵 \(\Phi\) 下，线性回归模型 (9.13) 的负对数似然可以写为 \[ - \log p(Y \mid X, \theta) = \frac{1}{2\sigma^2} (y - \Phi \theta)^{\top}(y - \Phi \theta) + \text{const}. \tag{9.18} \] 将 (9.18) 与 “无特征” 模型的负对数似然 (9.10b) 对比，我们可以立刻看到，只需要用 \(\Phi\) 替换 \(X\) 即可。由于 \(X\) 和 \(\Phi\) 都与我们希望优化的参数 \(\theta\) 无关，我们立刻得到非线性特征下的线性回归问题的最大似然估计： \[ \theta_{\mathrm{ML}} = (\Phi^{\top}\Phi)^{-1} \Phi^{\top} y . \tag{9.19} \]

备注：当我们在没有特征的情况下工作时，我们要求 \(X^{\top}X\) 可逆，这在 \(\mathrm{rk}(X) = D\)（即 \(X\) 的列向量线性无关）时成立。在 (9.19) 中，我们同样要求 \(\Phi^{\top}\Phi \in \mathbb{R}^{K \times K}\) 可逆。这当且仅当 \(\mathrm{rk}(\Phi) = K\) 时成立。

例 9.5（最大似然多项式拟合） 考虑图 9.4(a) 中的数据集。该数据集由 \(N=10\) 对 \((x_n, y_n)\) 组成，其中

\[ x_n \sim U[-5, 5], \quad y_n = -\sin(x_n/5) + \cos(x_n) + \epsilon \] s且 \(\epsilon \sim \mathcal{N}(0, 0.2^2)\)。

我们使用最大似然估计拟合一个 4 次多项式，即参数 \(\theta_{\text{ML}}\) 由公式 (9.19) 给出。最大似然估计在任意测试点 \(x^*\) 的预测值为\(\phi^\top(x^*) \theta_{\text{ML}}.\)结果如图 9.4(b) 所示。

个人注：例子里，真实的数据生成函数是： \[ y_n = -\sin(x_n/5) + \cos(x_n) + \epsilon \] 其中 \(\epsilon \sim N(0, 0.2^2)\) 是噪声。注意，这里的 真实函数是非多项式的，它是由正弦和余弦函数组成的。

但是在这个例子中，作者选择用 4次多项式 来拟合这些数据： \[ y \approx \theta_0 + \theta_1 x + \theta_2 x^2 + \theta_3 x^3 + \theta_4 x^4 \]

噪声方差的估计

到目前为止，我们一直假设噪声方差 σ² 是已知的。然而，我们同样可以利用最大似然估计的原理来得到噪声方差的最大似然估计量 \(σ²_{ML}\)。为此，我们遵循标准步骤：写出对数似然函数，对 \(σ² > 0\) 求导，使其等于 0，并解方程。

个人注：噪声方差的最大似然估计告诉你模型拟合后，数据中剩余的不确定性有多大。

噪声方差的最大似然估计（MLE for noise variance）在统计建模中有非常直观的意义，它告诉我们 模型预测值与实际观测值之间偏差的典型大小，也就是 数据中的随机波动有多大。换句话说，它量化了模型拟合“剩余误差”的平均能量。剩余误差就是数据本身的随机波动。

意义：

1. 平均残差平方：它就是残差平方的平均值，也就是模型拟合后剩余误差的能量。
2. 不确定性度量：在高斯假设下，方差越大，说明数据波动越大，模型预测的不确定性越高。
3. 概率解释：在高斯模型下，\(\hat{\sigma}^2\) 最大化了观测数据出现的概率（即似然函数最大化）。

对数似然函数为： \[ \begin{align} \log p(Y \mid X, \theta, \sigma^2) &= \sum_{n=1}^N \log \mathcal{N}\!\left(y_n \mid \phi^\top(x_n)\theta, \sigma^2 \right) \tag{9.20a} \\[6pt] &= \sum_{n=1}^N \left( -\tfrac{1}{2}\log(2\pi) -\tfrac{1}{2}\log\sigma^2 -\tfrac{1}{2\sigma^2}(y_n - \phi^\top(x_n)\theta)^2 \right) \tag{9.20b} \\[6pt] &= -\tfrac{N}{2}\log\sigma^2 -\tfrac{1}{2\sigma^2} \underbrace{\sum_{n=1}^N (y_n - \phi^\top(x_n)\theta)^2 \,}_{=:s} + \text{const.} \tag{9.20c} \end{align} \]

对数似然函数关于 σ² 的偏导数为： \[ \begin{align} \frac{\partial \log p(Y \mid X, \theta, \sigma^2)}{\partial \sigma^2} &= -\frac{N}{2\sigma^2} + \frac{1}{2\sigma^4}s = 0 \tag{9.21a} \\[6pt] \Longleftrightarrow \quad \frac{N}{2\sigma^2} &= \frac{s}{2\sigma^4} \tag{9.21b} \end{align} \]

因此我们得到： \[ \sigma^2_{ML} = \frac{s}{N} = \frac{1}{N}\sum_{n=1}^N \big(y_n - \phi^\top(x_n)\theta \big)^2 . \tag{9.22} \] 因此，噪声方差的最大似然估计就是在输入位置 \(x_n\) 上，噪声自由的函数值 \(\phi^\top(x_n)\theta\) 与相应带噪观测值 \(y_n\) 之间平方差的经验均值。

9.2.2 线性回归中的过拟合

我们刚刚讨论了如何使用最大似然估计来将线性模型（例如多项式）拟合到数据上。我们可以通过计算模型产生的误差/损失来评估模型的质量。一种方法是计算负对数似然（公式 9.10b），我们通过最小化它来确定最大似然估计器。另一种方法是注意到噪声参数 \(\sigma^2\) 不是一个自由的模型参数，因此我们可以忽略 \(1/\sigma^2\) 的缩放，从而得到平方误差损失函数 \[ \|y - \Phi \theta\|^2. \]

个人注：取负对数似然（negative log-likelihood, NLL）得到： \[ -\log p(y \mid \theta, \sigma^2) = \frac{N}{2} \log (2\pi \sigma^2) + \frac{1}{2\sigma^2} \sum_{n=1}^{N} (y_n - (\Phi \theta)_n)^2 \] 注意这里有两部分：

1. \(\frac{N}{2} \log (2\pi \sigma^2)\) —— 与 \(\theta\) 无关
2. \(\frac{1}{2\sigma^2} \sum (y_n - (\Phi \theta)_n)^2\) —— 与 \(\theta\) 有关，但前面有 \(1/(2\sigma^2)\) 的缩放

为什么可以忽略 \(1/\sigma^2\)

• 当我们只对 \(\theta\) 求最大似然估计时，\(\sigma^2\) 被认为是一个已知常数或独立处理的参数。
• 对于优化 \(\theta\) 来说，乘以一个正的常数 \(1/\sigma^2\) 并不会改变最小点的位置。
• 因此，为了简化计算，我们可以忽略 \(1/\sigma^2\) 的缩放因子，直接最小化平方误差：

\[ \sum_{n=1}^{N} (y_n - (\Phi \theta)_n)^2 = \|y - \Phi \theta\|^2 \]

这就是为什么在线性回归中，平方误差损失等价于最大似然估计的原因。

除了使用平方损失外，我们通常使用均方根误差（root mean square error）（RMSE）： \[ \sqrt{\frac{1}{N}\|y - \Phi\theta\|^2} \;=\; \sqrt{\frac{1}{N} \sum_{n=1}^N (y_n - \phi(x_n)^\top \theta)^2}, \tag{9.23} \] 它具有以下优点： (a) 允许我们比较不同规模数据集上的误差； (b) 与观测的函数值 \(y_n\) 具有相同的量纲和单位。

例如，如果我们拟合一个模型，将邮政编码（\(x\) 由经纬度表示）映射到房价（\(y\) 的单位是欧元），那么 RMSE 的单位也是欧元，而平方误差的单位则是欧元平方（EUR\(^2\)）。如果我们选择在目标函数中包含原始负对数似然（公式 9.10b）中的因子 \(\sigma^2\)，那么最终目标将不再有单位；在前述例子中，目标函数将不再是以欧元或欧元平方为单位。

在模型选择中（见第 8.6 节），我们可以使用 RMSE（或负对数似然）来确定多项式的最佳次数，即找到使目标函数最小化的多项式次数 \(M\)。由于多项式次数是自然数，我们可以通过穷举搜索来枚举所有（合理的）\(M\) 值。对于一个大小为 \(N\) 的训练集，测试 \(0 \leq M \leq N-1\) 就足够了。

当 \(M < N\) 时，最大似然估计是唯一的；而当 \(M \geq N\) 时，模型的参数数量多于数据点数量，我们就需要解一个欠定的线性方程组（此时 (9.19) 中的 \(\Phi^\top \Phi\) 也将不可逆），因此会有无穷多个可能的最大似然估计解。

F9.5

图 9.5 展示了利用最大似然方法在图 9.4(a) 中包含 \(N = 10\) 个观测点的数据集上得到的一系列多项式拟合结果。我们可以注意到，低阶多项式（例如常数函数 \(M=0\) 或一次函数 \(M=1\)）对数据的拟合效果很差，因此对真实的底层函数表现力不足。对于 \(M = 3, \ldots, 6\) 的情况，拟合结果看起来合理，能够平滑地插值通过这些数据点。当阶数进一步增大时，我们发现多项式对数据的拟合越来越好。在极端情况下，当 \(M = N - 1 = 9\) 时，拟合函数会精确通过每一个数据点。然而，这类高阶多项式会出现剧烈振荡，无法很好地表示生成数据的真实底层函数，因而会产生严重的过拟合问题。

请记住，我们的目标是通过对新的（未见过的）数据进行准确预测，从而实现良好的泛化性能。为了更直观地了解泛化性能与多项式阶数 \(M\) 之间的关系，我们考虑了一个独立的测试集，该测试集由 200 个数据点组成，数据生成方式与训练集完全相同。测试输入点选取为区间 \([-5, 5]\) 内的 200 个等间隔点。对于每一个 \(M\)，我们分别在训练集和测试集上计算 RMSE (公式 9.23)。

观察测试误差（它是评价对应多项式泛化能力的一个定性指标），我们发现其初始阶段会下降；参见图 9.6（橙色曲线）。对于四阶多项式，测试误差相对较小，并且在 \(M=5\) 之前基本保持稳定。然而，从 \(M=6\) 开始，测试误差显著上升，高阶多项式的泛化能力非常差。在这个例子中，这一点也可以从图 9.5 中的最大似然拟合曲线清楚地看出。需要注意的是，训练误差（图 9.6 中的蓝色曲线）在多项式阶数增加时从不会上升。在本例中，最好的泛化效果（测试误差最小的点）出现在 \(M=4\) 的多项式上。

9.2.3 最大后验估计

Maximum A Posteriori Estimation, MAP

我们刚刚看到，最大似然估计（MLE）容易产生过拟合。当出现过拟合时，参数值的幅度往往会变得相对较大（Bishop, 2006）。为了减轻参数值过大的影响，我们可以在参数上引入一个先验分布 \(p(\theta)\)。先验分布明确地编码了哪些参数值在观测数据之前是合理的。例如，对于单个参数 \(\theta\)，如果我们选择一个高斯先验 \(p(\theta) = \mathcal{N}(0, 1)\)，就表示我们期望 \(\theta\) 的取值落在区间 \([-2, 2]\) 内（均值的两个标准差范围内）。一旦数据集 \(X, Y\) 可用，我们不再只最大化似然函数，而是寻找能够最大化后验分布 \(p(\theta | X, Y)\) 的参数。这一过程称为 最大后验估计（MAP 估计）。根据贝叶斯定理（参见第 6.3 节），给定训练数据 \(X, Y\)，参数 \(\theta\) 的后验分布为

\[ p(\theta | X, Y) = \frac{p(Y | X, \theta)p(\theta)}{p(Y | X)} . \tag{9.24} \] 由于后验分布显式依赖于参数的先验 \(p(\theta)\)，因此先验会影响我们得到的后验最大化参数向量。后面我们会更清楚地看到这一点。能够最大化 (9.24) 式后验分布的参数向量 \(\theta_{\text{MAP}}\) 就是 MAP 估计。

求解 MAP 估计的步骤与最大似然估计类似。我们先对后验取对数，得到对数后验： \[ \log p(\theta | X, Y) = \log p(Y | X, \theta) + \log p(\theta) + \text{const}, \tag{9.25} \] 其中常数项不依赖于 \(\theta\)。我们可以看到，对数后验 (9.25) 是对数似然项 \(\log p(Y | X, \theta)\) 与对数先验 \(\log p(\theta)\) 的和，因此 MAP 估计实际上是 先验（观测数据之前对参数的合理性假设）与数据驱动的似然函数之间的折中。

个人注：所以说先验分布是一种正则化手段。

为了求解 MAP 估计 \(\theta_{\text{MAP}}\)，我们需要最小化关于 \(\theta\) 的负对数后验： \[ \theta_{\text{MAP}} \in \arg \min_{\theta} \big\{ -\log p(Y | X, \theta) - \log p(\theta) \big\}. \tag{9.26} \] 其梯度为： \[ - \frac{d \log p(\theta | X, Y)}{d\theta} = - \frac{d \log p(Y | X, \theta)}{d\theta} - \frac{d \log p(\theta)}{d\theta}, \tag{9.27} \] 其中右边第一个项就是 (9.11c) 中负对数似然的梯度。若在参数 \(\theta\) 上使用共轭的高斯先验

\[ p(\theta) = \mathcal{N}(0, b^2 I), \] 则在 (9.13) 的线性回归设定下，负对数后验为： \[ - \log p(\theta | X, Y) = \frac{1}{2\sigma^2} (y - \Phi \theta)^\top (y - \Phi \theta) + \frac{1}{2b^2} \theta^\top \theta + \text{const}. \tag{9.28} \] 其中第一项来自对数似然，第二项来自对数先验。其关于参数 \(\theta\) 的梯度为： \[ - \frac{d \log p(\theta | X, Y)}{d\theta} = \frac{1}{\sigma^2}(\theta^\top \Phi^\top \Phi - y^\top \Phi) + \frac{1}{b^2} \theta^\top. \tag{9.29} \]

我们将通过令该梯度为 \(0^\top\)，并解出 \(\theta_{\text{MAP}}\)，来找到最大后验估计 \(\theta_{\text{MAP}}\)。我们得到：

$$ \[\begin{align} & \quad \quad \quad \frac{1}{\sigma^2}(\theta^\top \Phi^\top \Phi - y^\top \Phi) + \frac{1}{b^2}\theta^\top = 0^\top \tag{9.30a} \\ & \Longleftrightarrow \quad \theta^\top \left( \frac{1}{\sigma^2}\Phi^\top \Phi + \frac{1}{b^2}I \right) - \frac{1}{\sigma^2}y^\top \Phi = 0^\top \tag{9.30b} \\ & \Longleftrightarrow \quad \theta^\top \left( \Phi^\top \Phi + \frac{\sigma^2}{b^2}I \right) = y^\top \Phi \tag{9.30c} \\ & \Longleftrightarrow \quad \theta^\top = y^\top \Phi \left( \Phi^\top \Phi + \frac{\sigma^2}{b^2}I \right)^{-1} \tag{9.30d} \end{align}\] $$

因此，MAP 估计为（通过对最后一个等式两边取转置）： \[ \theta_{\text{MAP}} = \left( \Phi^\top \Phi + \frac{\sigma^2}{b^2} I \right)^{-1} \Phi^\top y . \tag{9.31} \] 将 (9.31) 中的 MAP 估计与 (9.19) 中的最大似然估计进行比较，我们可以看到，两者解的唯一差别在于逆矩阵中额外的项 \(\tfrac{\sigma^2}{b^2} I\)。这个附加项保证了 \[ \Phi^\top \Phi + \frac{\sigma^2}{b^2} I \] 是对称的并且严格正定的（即其逆矩阵存在，并且 MAP 估计是该线性方程组的唯一解）。此外它反映了正则化项的影响。

例 9.6 （多项式回归的 MAP 估计） 在 9.2.1 节的多项式回归例子中，我们对参数 \(\theta\) 施加一个高斯先验 \(p(\theta) = \mathcal{N}(0, I)\)，并根据公式 (9.31) 来确定 MAP 估计。在图 9.7 中，我们展示了 6 次多项式（左图）和 8 次多项式（右图）的最大似然估计与 MAP 估计。对于低阶多项式，先验（正则化项）的作用不显著，但在高阶多项式时，它能使函数保持相对光滑。虽然 MAP 估计可以在一定程度上缓解过拟合，但它并不是该问题的通用解决方案，因此我们需要一种更为系统的方法来应对过拟合。

个人注：正则化是针对同一个模型，能让参数靠近原点，而不是改变参数的个数（模型）。

9.2.4 最大后验估计作为正则化

MAP Estimation as Regularizatio

个人注：感觉标题翻译为 “最大后验估计作和正则化” 更容易理解一点。

除了在参数 \(θ\) 上放置一个先验分布外，还可以通过正则化来惩罚参数的幅度，从而减轻过拟合的影响。在正则化的最小二乘法中，我们考虑如下的损失函数： \[ \|y - \Phi \theta\|^2 + \lambda \|\theta\|^2_2 , \tag{9.32} \] 并使其对 \(θ\) 最小化（参见第 8.2.3 节）。其中，第一个项是数据拟合项（也称为失配项 misfit term），它与负对数似然成正比（见公式 (9.10b)）。第二个项称为正则项 regularizer，正则化参数 \(\lambda \geq 0\) 控制了正则化的“严格程度”。

备注： 在 (9.32) 中，我们不一定要使用欧几里得范数 \(\|\cdot\|_2\)，而是可以选择任意的 p-范数 \(\|\cdot\|_p\)。在实践中，较小的 \(p\) 值会导致更稀疏的解。这里，“稀疏”意味着许多参数值满足 \(\theta_d = 0\)，这对于变量选择也很有用。当 \(p = 1\) 时，正则项被称为 LASSO（最小绝对收缩与选择算子）(least absolute shrinkage and selection operator)，由 Tibshirani (1996) 提出。

在式 (9.32) 中的正则项 \(\lambda \lVert \theta \rVert_2^2\) 可以解释为一个高斯先验的负对数形式，我们在 MAP 估计中会使用它；参见 (9.26)。更具体地说，若我们采用高斯先验\(p(\theta) = \mathcal{N}(0, b^2 I),\)那么对应的负对数高斯先验为 \[ - \log p(\theta) = \frac{1}{2b^2} \lVert \theta \rVert_2^2 + \text{const}, \tag{9.33} \] 因此，当 \(\lambda = \tfrac{1}{2b^2}\) 时，正则项与负对数高斯先验是完全相同的。考虑到在 (9.32) 中的正则化最小二乘损失函数由与负对数似然加负对数先验密切相关的项组成，因此，当我们最小化这个损失函数时，得到的解与 (9.31) 中的 MAP 估计非常相似就不足为奇了。更具体地，最小化正则化最小二乘损失函数得到

\[ \theta_{\text{RLS}} = (\Phi^\top \Phi + \lambda I)^{-1} \Phi^\top y, \tag{9.34} \] 当 \(\lambda = \tfrac{\sigma^2}{b^2}\) 时，这个解与 (9.31) 中的 MAP 估计完全一致，其中 \(\sigma^2\) 是噪声方差，而 \(b^2\) 是各向同性高斯先验\(p(\theta) = \mathcal{N}(0, b^2 I)\)的方差。

到目前为止，我们已经介绍了使用最大似然估计（MLE）和最大后验估计（MAP）的参数估计方法，在这些方法中，我们找到的是使目标函数（似然函数或后验概率）最大化的点估计 \(\theta^*\)。我们看到，无论是最大似然估计还是 MAP 估计，都可能导致过拟合。接下来的部分，我们将讨论贝叶斯线性回归，在这种方法中，我们使用贝叶斯推断（见第 8.4 节）来得到未知参数的后验分布，并进一步用它来进行预测。更具体地说，在预测时，我们会在所有合理的参数集合上进行平均，而不是仅仅依赖某一个点估计。

9.3 贝叶斯线性回归

在之前的内容中，我们研究了线性回归模型，并通过最大似然或 MAP 方法来估计模型参数 \(\theta\)。我们发现，MLE 在小数据量的情况下可能会导致严重的过拟合。MAP 方法通过在参数上引入一个先验分布，起到了正则化的作用，从而在一定程度上缓解了这一问题。

贝叶斯线性回归将“参数先验”的思想进一步推进：它甚至不再尝试计算参数的点估计，而是通过参数的完整后验分布来进行预测。换句话说，我们并不拟合某一个具体的参数，而是根据后验分布，对所有可能的参数设定进行加权平均。

9.3.1 模型

在贝叶斯线性回归中，我们考虑如下模型： \[ \begin{align} & \text{prior } \;\quad \quad p(\theta) = \mathcal{N}(m_0, S_0), \\ & \text{likelihood } \; p(y \mid x, \theta) = \mathcal{N}\big(y \mid \phi^{\top}(x)\theta, \sigma^2\big). \tag{9.35} \end{align} \] 这里我们显式地对参数向量 \(\theta\) 施加一个高斯先验 \(p(\theta) = \mathcal{N}(m_0, S_0)\)，从而把参数向量转变为一个随机变量。

这使我们能够写出对应的图模型，如图 9.8 所示，在该图中我们将 \(\theta\) 的高斯先验的参数明确表示出来。完整的概率模型（即观测变量与未观测变量的联合分布，分别是 \(y\) 和 \(\theta\)）可以写为： \[ p(y, \theta \mid x) = p(y \mid x, \theta)\, p(\theta). \tag{9.36} \]

9.3.2 先验预测

在实践中，我们通常并不是特别关心参数值 \(\theta\) 本身。相反，我们更关注的是利用这些参数值所做出的预测。在贝叶斯设定下，当我们进行预测时，会基于参数分布，对所有可能的参数取平均。更具体地说，为了在输入 \(x_*\) 处进行预测，我们需要对 \(\theta\) 积分，从而得到：

\[ p(y_* \mid x_*) = \int p(y_* \mid x_*, \theta)p(\theta)\, d\theta = \mathbb{E}_\theta\!\left[p(y_* \mid x_*, \theta)\right], \tag{9.37} \] 这可以解释为：对所有依据先验分布 \(p(\theta)\) 而合理的参数 \(\theta\)，其条件预测 \(y_* \mid x_*, \theta\) 的平均。注意，使用先验分布进行预测时，只需要给定输入 \(x_*\)，而不需要任何训练数据。

在我们的模型 (9.35) 中，我们为 \(\theta\) 选择了一个共轭（高斯）先验，因此预测分布同样是高斯分布（并且可以解析地计算）。具体来说，若先验分布为\(p(\theta) = \mathcal{N}(m_0, S_0),\)则预测分布为 \[ p(y_* \mid x_*) = \mathcal{N}\!\Big(\phi^\top(x_*) m_0, \, \phi^\top(x_*) S_0 \phi(x_*) + \sigma^2\Big), \tag{9.38} \] 这里我们利用了以下事实：

1. 由于共轭性（见第 6.6 节）和高斯分布的边缘化性质（见第 6.5 节），预测分布仍然是高斯分布；

2. 高斯噪声独立，因此 \[ \mathrm{Var}[y_*] = \mathrm{Var}_\theta[\phi^\top(x_*)\theta] + \mathrm{Var}_\epsilon[\epsilon], \tag{9.39} \]

3. \(y_*\) 是 \(\theta\) 的线性变换，因此可以直接应用 (6.50) 和 (6.51) 中的规则，解析计算预测的均值和协方差。

在公式 (9.38) 中，预测方差中的项 \(\phi^\top(x_*) S_0 \phi(x_*)\) 显式体现了与参数 \(\theta\) 相关的不确定性，而 \(\sigma^2\) 则反映了测量噪声带来的不确定性。如果我们更关心的是预测无噪声的函数值\(f(x_*) = \phi^\top(x_*)\theta,\)而不是被噪声污染的目标 \(y_*\)，则我们得到： \[ p(f(x_*)) = \mathcal{N}\!\Big(\phi^\top(x_*) m_0, \, \phi^\top(x_*) S_0 \phi(x_*)\Big), \tag{9.40} \] 与 (9.38) 的唯一区别在于预测方差中去掉了噪声方差 \(\sigma^2\)。

备注（函数分布）： 由于我们可以用一组样本 \(\{\theta_i\}\) 来表示参数分布 \(p(\theta)\)，并且每个样本 \(\theta_i\) 都对应着一个函数 \[ f_i(\cdot) = \theta_i^\top \phi(\cdot)， \] 因此参数分布 \(p(\theta)\) 诱导出了一个关于函数的分布 \(p(f(\cdot))\)。这里符号 \((\cdot)\) 明确表示了函数关系。

例 9.7 （函数上的先验）

考虑一个五次多项式的贝叶斯线性回归问题。我们选择参数先验为\(p(\theta) = \mathcal{N}(0, \tfrac{1}{4}I)。\)

图 9.9 展示了该参数先验所诱导的函数先验分布（阴影区域：深灰色表示 67% 置信区间；浅灰色表示 95% 置信区间），并给出了一些从该先验中采样得到的函数。

一个函数样本的生成方式是：首先从先验中采样一个参数向量 \[ \theta_i \sim p(\theta)， \] 然后计算 \[ f_i(\cdot) = \theta_i^\top \phi(\cdot)。 \] 我们使用了 200 个输入点 \(x_* \in [-5, 5]\)，并将其映射到特征函数 \(\phi(\cdot)\)。图 9.9 中的不确定性（阴影区域）完全来自参数的不确定性，因为这里我们考虑的是无噪声的预测分布 (9.40)。

到目前为止，我们一直在基于参数先验 \(p(\theta)\) 来计算预测。然而，当我们拥有参数后验（给定一些训练数据 \(X, Y\)）时，预测与推断的原理与公式 (9.37) 相同——只是需要将先验 \(p(\theta)\) 替换为后验 \(p(\theta \mid X,Y)\)。在接下来的内容中，我们将详细推导后验分布，然后再利用它进行预测。

9.3.3 后验分布

给定一组训练输入 \(\{x_n \in \mathbb{R}^D\}\) 及其对应的观测值 \(\{y_n \in \mathbb{R}\}, n=1,\dots,N\)，我们使用贝叶斯定理计算参数的后验分布： \[ p(\theta \mid X,Y) = \frac{p(Y \mid X,\theta) p(\theta)}{p(Y \mid X)}, \tag{9.41} \] 其中：\(X\) 是训练输入的集合； \(Y\) 是对应训练目标的集合；\(p(Y \mid X,\theta)\) 是似然函数；\(p(\theta)\) 是参数先验；

\[ p(Y \mid X) = \int p(Y \mid X,\theta)p(\theta)\, d\theta = \mathbb{E}_\theta[p(Y \mid X,\theta)] \tag{9.42} \]

称为 边际似然（证据）(marginal likelihood/evidence)，它与参数 \(\theta\) 无关，并确保后验归一化（积分为 1）。我们可以把边际似然看作是所有可能参数设置下（依据先验分布加权）的似然的平均。

定理 9.1 （参数后验） 在我们的模型 (9.35) 中，参数后验 (9.41) 可以解析地计算为：

\[ \begin{align} p(\theta \mid X,Y) & = \mathcal{N}(\theta \mid m_N, S_N), \tag{9.43a} \\ S_N & = (S_0^{-1} + \sigma^{-2}\Phi^\top \Phi)^{-1}, \tag{9.43b}\\ m_N & = S_N \big(S_0^{-1} m_0 + \sigma^{-2}\Phi^\top y \big), \tag{9.43c} \end{align} \]

下标 \(N\) 表示训练集大小。

证明

贝叶斯定理告诉我们，后验分布 \(p(\theta \mid X,Y)\) 正比于似然 \(p(Y \mid X,\theta)\) 与先验 \(p(\theta)\) 的乘积： \[ \begin{align} & \text{后验: }\qquad p(\theta \mid X,Y) = \frac{p(Y \mid X,\theta)p(\theta)}{p(Y \mid X)} \tag{9.44a} \\ & \text{似然: } \qquad p(Y \mid X,\theta) = \mathcal{N}(y \mid \Phi\theta, \sigma^2 I) \tag{9.44b} \\ & \text{先验: } \qquad p(\theta) = \mathcal{N}(\theta \mid m_0, S_0) \tag{9.44c} \end{align} \] 与其直接研究先验与似然的乘积，我们可以将问题转化到对数空间，并通过“配方”来求解后验的均值和协方差。先验对数与似然对数之和为：

\[ \begin{align} & \log \mathcal{N}(y \mid \Phi\theta, \sigma^2 I) + \log \mathcal{N}(\theta \mid m_0, S_0) \tag{9.45a} \\ & = -\tfrac{1}{2}\sigma^{-2}(y-\Phi\theta)^\top(y-\Phi\theta) \;-\; \tfrac{1}{2}(\theta - m_0)^\top S_0^{-1} (\theta - m_0) + \text{const}, \tag{9.45b} \end{align} \]

其中常数项与 \(\theta\) 无关，我们在后续忽略它。继续展开 (9.45b)，得到： \[ -\tfrac{1}{2}\sigma^{-2} y^\top y \;+\; \sigma^{-2} y^\top \Phi \theta \;-\; \tfrac{1}{2}\theta^\top \sigma^{-2}\Phi^\top \Phi \theta \;-\; \tfrac{1}{2}\theta^\top S_0^{-1}\theta \;+\; m_0^\top S_0^{-1}\theta \;-\; \tfrac{1}{2}m_0^\top S_0^{-1} m_0 + \text{const} \tag{9.46a} \] 整理后得到： \[ = -\tfrac{1}{2}\theta^\top(\sigma^{-2}\Phi^\top \Phi + S_0^{-1})\theta \;+\; (\sigma^{-2}\Phi^\top y + S_0^{-1}m_0)^\top \theta + \text{const}, \tag{9.46b} \] 显然这是关于 \(\theta\) 的二次型。未归一化的对数后验是一个负二次型，因此后验必然是高斯分布：

\[ \begin{align} & p(\theta \mid X,Y) = \exp(\log p(\theta \mid X,Y)) \propto \exp(\log p(Y \mid X,\theta) + \log p(\theta)) \tag{9.47a} \\ & \propto \exp\!\left(-\tfrac{1}{2}\theta^\top (\sigma^{-2}\Phi^\top \Phi + S_0^{-1})\theta \;+\; (\sigma^{-2}\Phi^\top y + S_0^{-1}m_0)^\top \theta \right), \tag{9.47b} \end{align} \]

由此可见，后验分布为高斯分布，其均值与协方差正是 (9.43b) 和 (9.43c) 中所给。剩下的任务就是将这个（未标准化的）高斯分布写成与\(\mathcal{N}(\theta \mid m_N, S_N)\)成比例的形式，也就是说，我们需要确定均值 \(m_N\) 和协方差矩阵 \(S_N\)。为此，我们使用配方（completing the squares）的技巧。所需的对数后验分布为

\[ \begin{align} \log \mathcal{N}(\theta \mid m_N, S_N) & = -\tfrac{1}{2} (\theta - m_N)^\top S_N^{-1} (\theta - m_N) + \text{const} \tag{9.48a} \\ & = -\tfrac{1}{2}\theta^\top S_N^{-1}\theta - m_N^\top S_N^{-1}\theta + \tfrac{1}{2} m_N^\top S_N^{-1} m_N. \tag{9.48b} \end{align} \]

在这里，我们将二次型 \((\theta - m_N)^\top S_N^{-1} (\theta - m_N)\) 拆分成三类项：关于 \(\theta\) 的二次项（蓝色），关于 \(\theta\) 的一次项（橙色），常数项（黑色）。(个人注：原书的公式是标了颜色。)这样我们就可以通过在 (9.46b) 和 (9.48b) 中对比相同颜色的项来得到 \(S_N\) 和 \(m_N\)：

\[ \begin{align} & \quad S_N^{-1} = \Phi^\top \sigma^{-2} I \Phi + S_0^{-1} \tag{9.49a} \\ \Longleftrightarrow & \quad S_N = \left(\sigma^{-2}\Phi^\top \Phi + S_0^{-1}\right)^{-1} \tag{9.49b} \end{align} \]

并且有

\[ \begin{align} & \quad m_N^\top S_N^{-1} = \left(\sigma^{-2}\Phi^\top y + S_0^{-1} m_0 \right)^\top \tag{9.50a} \\ \Longleftrightarrow \quad & \quad m_N = S_N \left(\sigma^{-2}\Phi^\top y + S_0^{-1} m_0 \right). \tag{9.50b} \end{align} \]

注释（配方法的一般方法）。如果我们有如下形式的等式： \[ x^\top A x - 2 a^\top x + \text{const}_1 , \tag{9.51} \]

其中 \(A\) 是对称且正定的，我们希望将其改写成如下形式：

\[ (x - \mu)^\top \Sigma (x - \mu) + \text{const}_2 , \tag{9.52} \]

那么我们可以通过以下设定来实现：

\[ \Sigma := A, \tag{9.53} \]

\[ \mu := \Sigma^{-1} a, \tag{9.54} \]

并且

\[ \text{const}_2 = \text{const}_1 - \mu^\top \Sigma \mu. \quad \]

我们可以看到，(9.47b) 中指数里的项正是 (9.51) 的形式，其中：

\[ \begin{align} A & := \sigma^{-2}\Phi^\top \Phi + S_0^{-1}, \tag{9.55}\\ a & := \sigma^{-2}\Phi^\top y + S_0^{-1} m_0. \tag{9.56} \end{align} \]

由于在像 (9.46a) 这样的等式中，识别 \(A, a\) 可能比较困难，因此把这些等式改写为 (9.51) 的形式通常是有帮助的。这样可以把二次项、一次项和常数项分离开，从而简化所需解的推导。

9.3.4 后验预测

在式 (9.37) 中，我们利用参数先验 \(p(\theta)\) 计算了测试输入 \(x_*\) 处的预测分布 \(y_*\)。原则上，使用参数后验 \(p(\theta \mid X, Y)\) 进行预测并没有本质上的区别，因为在我们的共轭模型中，先验和后验都是高斯分布（只是参数不同）。因此，沿用第 9.3.2 节中的相同推理，我们得到（后验）预测分布：

\[ \begin{align} p(y_* \mid X, Y, x_*) & = \int p(y_* \mid x_*, \theta)p(\theta \mid X, Y)\, d\theta \tag{9.57a} \\ & = \int \mathcal{N}\!\big(y_* \mid \phi^\top(x_*) \theta, \sigma^2\big)\, \mathcal{N}\!\big(\theta \mid m_N, S_N\big)\, d\theta \tag{9.57b} \\ & = \mathcal{N}\!\Big(y_* \mid \phi^\top(x_*) m_N,\; \phi^\top(x_*) S_N \phi(x_*) + \sigma^2\Big). \tag{9.57c} \end{align} \]

其中，项 \(\phi^\top(x_*) S_N \phi(x_*)\) 反映了与参数 \(\theta\) 相关的后验不确定性。注意，\(S_N\) 依赖于训练输入（通过 \(\Phi\)）；参见式 (9.43b)。预测均值 \(\phi^\top(x_*) m_N\) 与最大后验估计（MAP 估计）\(\theta_{\text{MAP}}\) 的预测结果一致。

注释（边际似然与后验预测分布）。通过替换式 (9.57a) 中的积分，预测分布也可以等价地写为： \[ \mathbb{E}_{\theta \mid X,Y}[p(y_* \mid x_*, \theta)], \]

其中期望是关于参数后验 \(p(\theta \mid X,Y)\) 进行的。这样表述后验预测分布凸显了其与边际似然 (9.42) 的紧密联系。两者的关键区别在于：

1. 边际似然可以被看作是在预测训练目标 \(y\)，而不是测试目标 \(y_*\)；
2. 边际似然是对参数先验进行平均，而后验预测分布是对参数后验进行平均。

注释（无噪声函数值的均值与方差）。在很多情况下，我们并不关心（有噪声的）观测 \(y_*\) 的预测分布 \(p(y_* \mid X,Y,x_*)\)，而是希望得到（无噪声的）函数值 \[ f(x_*) = \phi^\top(x_*) \theta \]

的分布。我们可以利用均值与方差的性质计算其对应矩，从而得到：

\[ \begin{align} \mathbb{E}[f(x_*) \mid X,Y] & = \mathbb{E}_\theta[\phi^\top(x_*) \theta \mid X,Y] = \phi^\top(x_*) \mathbb{E}_\theta[\theta \mid X,Y] \tag{9.58} \\ & = \phi^\top(x_*) m_N = m_N^\top \phi(x_*), \end{align} \]

\[ \begin{align} \mathbb{V}_\theta[f(x_*) \mid X,Y] & = \mathbb{V}_\theta[\phi^\top(x_*) \theta \mid X,Y] \tag{9.59} \\ & = \phi^\top(x_*) \mathbb{V}_\theta[\theta \mid X,Y] \phi(x_*) \\ & = \phi^\top(x_*) S_N \phi(x_*). \end{align} \]

我们看到，由于噪声的均值为 0，预测均值与含噪观测的预测均值相同；预测方差仅在于多了一个 σ² 的差异，这个 σ² 就是测量噪声的方差：当我们预测含噪的函数值时，需要将 σ² 作为不确定性的来源纳入考虑，但在无噪声预测中则不需要这一项。在这里，唯一剩下的不确定性来自于参数的后验分布。

备注（函数上的分布）：我们将参数 θ 积分消去的事实，诱导出了函数上的一个分布：如果我们从参数后验分布 \(p(\theta \mid X,Y)\) 中采样 \(\theta_i\)，那么我们得到的就是一个函数实现 \(\theta_i^\top \phi(\cdot)\)。这个函数分布的均值函数（即所有期望函数值 \(\mathbb{E}_\theta[f(\cdot)\mid \theta,X,Y]\) 的集合）为 \(m_N^\top \phi(\cdot)\)。其（边际）方差，即函数 \(f(\cdot)\) 的方差，则由 \(\phi^\top(\cdot) S_N \phi(\cdot)\) 给出。

例 9.8（函数上的后验分布） 让我们重新考察多项式次数为 5 的贝叶斯线性回归问题。我们选择参数先验为

\[ p(\theta) = \mathcal{N}(0, \tfrac{1}{4}I). \]

图 9.9 展示了由该参数先验诱导的函数先验分布，以及从该先验中采样得到的一些函数。

图 9.10 展示了通过 贝叶斯线性回归 得到的函数后验分布。面板 (a) 显示了训练数据集；面板 (b) 展示了函数的后验分布，其中包括通过最大似然估计 (MLE) 和最大后验估计 (MAP) 得到的函数。在贝叶斯线性回归中，通过 MAP 得到的函数也对应于 后验均值函数。面板 (c) 展示了在该函数后验分布下可能出现的一些函数实例（样本）。

图 9.11 展示了由参数后验分布诱导出的若干函数后验分布。对于不同的多项式阶数 \(M\)，左侧面板展示了最大似然函数 \(\theta^\top_{\text{ML}}\phi(\cdot)\)、MAP 函数 \(\theta^\top_{\text{MAP}}\phi(\cdot)\)（它与后验均值函数相同），以及通过贝叶斯线性回归得到的 67% 和 95% 的预测置信区间（用阴影区域表示）。

右侧面板展示了函数后验分布中的一些采样：在这里，我们从参数后验分布中采样参数 \(\theta_i\)，并计算函数 \(\phi^\top(x^*)\theta_i\)，这对应于后验分布下的某个函数实例。对于低阶多项式，参数后验分布限制了参数的变化范围，因此采样得到的函数几乎相同。当我们通过增加更多参数使模型更灵活（即更高阶的多项式）时，这些参数在后验分布中受到的约束不足，因而采样的函数可以被轻易区分开来。我们还可以在左侧对应的面板中看到，不确定性随之增加，特别是在边界区域。

尽管对于七阶多项式，MAP 估计能够给出一个合理的拟合，但贝叶斯线性回归模型还额外告诉我们：后验不确定性非常大。在将这些预测用于决策系统时，这一点信息可能至关重要，因为错误的决策可能会带来严重后果（例如在强化学习或机器人应用中）。

个人注：贝叶斯线性回归还给出不确定性的情况！

9.3.5 计算边际似然

在 第 8.6.2 节 中，我们强调了边际似然在贝叶斯模型选择中的重要性。接下来，我们将在带有共轭高斯先验的参数化贝叶斯线性回归模型中计算边际似然，也就是本章所讨论的设定。

回顾一下，我们考虑如下的生成过程： \[ \begin{align} \theta & \sim \mathcal{N}(m_0, S_0) \tag{9.60a} \\ y_n | x_n, \theta & \sim \mathcal{N}(x_n^\top \theta, \sigma^2), \quad n = 1, \ldots, N. \tag{9.60b} \end{align} \]

边际似然由下式给出：

\[ \begin{align} p(Y|X) & = \int p(Y|X, \theta)\,p(\theta)\,d\theta \tag{9.61a} \\ & = \int \mathcal{N}(y|X\theta, \sigma^2 I)\,\mathcal{N}(\theta|m_0, S_0)\, d\theta, \tag{9.61b} \end{align} \]

其中我们对模型参数 \(\theta\) 进行了积分。我们分两步计算边际似然：

1. 边际似然是高斯分布

在 第 6.5.2 节 中，我们知道：

1. 两个高斯随机变量的乘积仍然是（未归一化的）高斯分布；
2. 高斯随机变量的线性变换仍然是高斯分布。

在 (9.61b) 中，我们需要将 \(\mathcal{N}(y|X\theta, \sigma^2 I)\) 转换为某个形式 \(\mathcal{N}(\theta|\mu, \Sigma)\)。一旦完成转换，该积分就可以闭式解出，结果就是这两个高斯分布乘积的归一化常数。而归一化常数本身呈现高斯分布形式（见公式 (6.76)）。

2. 均值和协方差

我们利用仿射变换下随机变量的均值和协方差的标准结果（见第 6.4.4 节）来计算边际似然的均值与协方差矩阵。

边际似然的均值为： \[ \mathbb{E}[Y|X] = \mathbb{E}_{\theta,\epsilon}[X\theta + \epsilon] = X\mathbb{E}_\theta[\theta] = Xm_0. \tag{9.62} \]

其中 \(\epsilon \sim \mathcal{N}(0, \sigma^2 I)\) 是一个独立同分布的随机向量。

协方差矩阵为： \[ \begin{align} \mathrm{Cov}[Y|X] & = \mathrm{Cov}_{\theta,\epsilon}[X\theta + \epsilon] = \mathrm{Cov}_\theta[X\theta] + \sigma^2 I \tag{9.63a}\\ & = X\mathrm{Cov}_\theta[\theta]X^\top + \sigma^2 I = XS_0X^\top + \sigma^2 I. \tag{9.63b} \end{align} \]

因此，边际似然为：

\[ \begin{align} p(Y|X) & = (2\pi)^{-\frac{N}{2}} \det(XS_0X^\top + \sigma^2 I)^{-\frac{1}{2}} \exp\!\Bigg(-\tfrac{1}{2}(y - Xm_0)^\top (XS_0X^\top + \sigma^2 I)^{-1} (y - Xm_0)\Bigg) \tag{9.64a} \\ & = \mathcal{N}\big(y \,|\, Xm_0,\, XS_0X^\top + \sigma^2 I\big) \tag{9.64b} \end{align} \]

鉴于其与 后验预测分布 的紧密联系（参见本节前面的备注 “边际似然与后验预测分布”），边际似然的函数形式并不令人意外。

9.4 最大似然作为正交投影

在通过大量代数运算推导出最大似然和 MAP（最大后验）估计之后，我们现在给出最大似然估计的一种几何解释。考虑一个简单的线性回归情形：

\[ y = x\theta + \epsilon,\quad \epsilon \sim \mathcal{N}(0, \sigma^2), \tag{9.65} \]

其中我们考虑通过原点的线性函数 \(f: \mathbb{R} \to \mathbb{R}\)（为了简洁这里省略了特征）。参数 \(\theta\) 决定了直线的斜率。图 9.12(a) 显示了一个一维数据集。

给定训练数据集 \(\{(x_1,y_1),\dots,(x_N,y_N)\}\)，回顾 9.2.1 节的结果，可以得到斜率参数的最大似然估计：

\[ \theta_{\text{ML}} = (X^\top X)^{-1}X^\top y = \frac{X^\top y}{X^\top X} \in \mathbb{R}, \tag{9.66} \]

其中 \[ X = [x_1,\dots,x_N]^\top \in \mathbb{R}^N, \; y = [y_1,\dots,y_N]^\top \in \mathbb{R}^N \]

个人注：相当于在一个\(n\)维空间中，向量\(y\)投影到向量\(X\)上，因为投影\(\bar y\)是落在向量\(X\)上，故\(\bar{y}_i= ax_i\)；从而拟合的曲线在二维空间中是\(\bar y=ax\)；具体详见 “3.8.1 一维子空间（线）上的投影” 的公式（ 3.42）

这意味着对于训练输入 \(X\)，我们得到对训练目标的最优（最大似然）重建为： \[ X\theta_{\text{ML}} = X\frac{X^\top y}{X^\top X} = \frac{XX^\top}{X^\top X}y, \tag{9.67} \]

即得到在 \(y\) 和 \(X\theta\) 之间最小化平方误差的近似。

由于我们正在寻找方程 \(y = X\theta\) 的解，可以将线性回归看作是解线性方程组的问题。因此，我们可以联系第 2、3 章中讨论的线性代数和解析几何概念。特别是仔细观察 (9.67)，可以看到在 (9.65) 的例子中，最大似然估计 \(\theta_{\text{ML}}\) 实际上是将 \(y\) 正交投影到由 \(X\) 张成的一维子空间上。回顾 3.8 节关于正交投影的结果，我们可以识别出 \(\tfrac{XX^\top}{X^\top X}\) 是投影矩阵，\(\theta_{\text{ML}}\) 是投影到 \(X\) 张成的一维子空间上的坐标，而 \(X\theta_{\text{ML}}\) 则是 \(y\) 在该子空间上的正交投影。

因此，最大似然解也提供了几何上最优的解，它找到由 \(X\) 张成的子空间中“最接近”观测值 \(y\) 的向量，其中“最接近”指的是函数值 \(y_n\) 与 \(x_n\theta\) 之间平方距离最小。这一过程是通过正交投影实现的。图 9.12(b) 显示了带噪观测值投影到子空间上的情况，它最小化了原始数据集与投影之间的平方距离（注意 \(x\)-坐标是固定的），这正对应于最大似然解。

在一般的线性回归情形下：

\[ y = \phi^\top(x)\theta + \epsilon,\quad \epsilon \sim \mathcal{N}(0,\sigma^2), \tag{9.68} \]

其中 \(\phi(x) \in \mathbb{R}^K\) 是向量值特征，我们同样可以将最大似然结果解释为：

\[ y \approx \Phi \theta_{\text{ML}}, \tag{9.69} \]

\[ \theta_{\text{ML}} = (\Phi^\top \Phi)^{-1}\Phi^\top y, \tag{9.70} \]

即投影到由特征矩阵 \(\Phi\) 的列所张成的 \(\mathbb{R}^N\) 的 \(K\) 维子空间上；参见 3.8.2 节。

个人注：特征矩阵详见 “9.2.1 极大似然估计”的公式（9.16）。

如果我们用来构造特征矩阵 \(\Phi\) 的特征函数 \(\phi_k\) 是正交归一的（见 3.7 节），那么我们得到一种特殊情况，此时 \(\Phi\) 的列构成一个正交归一基（见 3.5 节），因而：

\[ \Phi^\top \Phi = I. \]

这将导致投影：

\[ \Phi(\Phi^\top \Phi)^{-1}\Phi^\top y = \Phi \Phi^\top y = \left ( \sum_{k=1}^K \phi_k \phi_k^\top \right ) y, \tag{9.71} \]

因此最大似然投影就只是将 \(y\) 投影到各个基向量 \(\phi_k\)（即 \(\Phi\) 的列）上的投影之和。进一步地，由于基是正交的，不同特征之间的耦合消失了。信号处理中许多常见的基函数（如小波和傅里叶基）就是正交基函数。

当基不是正交的时，可以通过 Gram-Schmidt 过程将一组线性无关的基函数转换为正交基；参见 3.8.3 节和 Strang (2003)。

9.5 延伸阅读

在本章中，我们讨论了针对高斯似然以及模型参数的共轭高斯先验的线性回归。这使得贝叶斯推断可以得到闭式解。然而，在一些应用中，我们可能希望选择不同的似然函数。例如，在二分类问题中，我们只观察到两个可能的（类别型）结果，此时使用高斯似然并不合适。相反，我们可以选择 伯努利（Bernoulli）似然，它会返回预测标签为 1（或 0）的概率。关于分类问题的深入介绍，可参见 Barber (2012)、Bishop (2006) 和 Murphy (2012) 的著作。另一个非高斯似然重要的例子是 计数数据。计数是非负整数，在这种情况下，二项分布或泊松分布的似然比高斯似然更合适。所有这些例子都属于 广义线性模型（GLM） 的范畴，GLM 是线性回归的一种灵活推广，允许响应变量具有非高斯的误差分布。GLM 推广了线性回归，使得线性模型可以通过一个光滑且可逆的函数 \(\sigma(\cdot)\)（可能是非线性的）与观测值相关联，因此有

\[ y = \sigma(f(x)), \quad f(x) = \theta^\top \phi(x), \]

其中 \(f(x)\) 就是 (9.13) 中的线性回归模型。换句话说，可以将广义线性模型看作函数复合：

\[ y = \sigma \circ f, \]

其中 \(f\) 是线性回归模型，\(\sigma\) 是激活函数（activation function）。需要注意的是，尽管我们称其为“广义线性模型”，但输出 \(y\) 不再是参数 \(\theta\) 的线性函数。

在逻辑回归中，我们选择逻辑 sigmoid 函数

\[ \sigma(f) = \frac{1}{1 + \exp(-f)} \in [0,1], \]

它可以解释为伯努利随机变量 \(y \in \{0,1\}\) 中观测到 \(y=1\) 的概率。函数 \(\sigma(\cdot)\) 被称为 传递函数（transfer function）或激活函数（activation function），其逆函数被称为 规范链接函数（canonical link function）。从这个角度也可以清楚地看到，广义线性模型是（深度）前馈神经网络的构建模块：如果我们考虑一个广义线性模型

\[ y = \sigma(Ax + b), \]

其中 \(A\) 是权重矩阵，\(b\) 是偏置向量，我们就可以将该广义线性模型识别为一个带有激活函数 \(\sigma(\cdot)\) 的单层神经网络。我们可以递归地复合这些函数：

\[ x_{k+1} = f_k(x_k), \quad f_k(x_k) = \sigma_k(A_k x_k + b_k), \quad k = 0, \dots, K-1, \tag{9.72} \]

其中 \(x_0\) 是输入特征，\(x_K = y\) 是观测输出，这样 \(f_{K-1}\circ \cdots \circ f_0\) 就是一个 \(K\) 层深度神经网络。因此，这个深度神经网络的基本构建模块正是 (9.72) 中定义的广义线性模型。神经网络（Bishop, 1995; Goodfellow 等, 2016）比线性回归模型更具表现力和灵活性。然而，最大似然参数估计是一个 非凸优化问题，在完全贝叶斯的框架下对参数进行边缘化在解析上是不可行的。

我们曾简要提到过，参数上的分布会诱导出回归函数上的分布。高斯过程（Gaussian process）（Rasmussen 和 Williams, 2006）是一种以函数分布为核心的回归模型。高斯过程不是在参数上放置分布，而是直接在函数空间上放置分布，从而避免了“绕道”参数空间。为了实现这一点，高斯过程利用了 核技巧（kernel trick）（Schölkopf 和 Smola, 2002），它允许我们仅通过查看输入 \(x_i, x_j\)，就能计算函数值 \(f(x_i), f(x_j)\) 之间的内积。高斯过程与贝叶斯线性回归和支持向量回归密切相关，但也可以解释为一个单隐层的贝叶斯神经网络，其单元数趋于无穷（Neal, 1996; Williams, 1997）。优秀的高斯过程入门资料可参见 MacKay (1998) 和 Rasmussen 与 Williams (2006)。

在本章的讨论中，我们专注于高斯先验，因为它能使线性回归模型的推断得到闭式解。然而，即使在具有高斯似然的回归设定中，我们也可能选择非高斯先验。考虑一种情形：输入为 \(x \in \mathbb{R}^D\)，但训练集规模很小，\(N \ll D\)。这意味着回归问题是 欠定的。在这种情况下，我们可以选择一种强制稀疏性的参数先验，即尽可能将更多参数压缩为 0 的先验（变量选择）。这种先验比高斯先验提供了更强的正则化，通常会带来更高的预测准确率和更好的模型可解释性。

拉普拉斯先验（Laplace prior） 是一个常用的例子。带有拉普拉斯先验的线性回归模型等价于带 \(L_1\) 正则化的线性回归（LASSO）（Tibshirani, 1996）。拉普拉斯分布在零点处尖锐（其一阶导数不连续），并且将概率质量更集中在零附近，相比高斯分布更鼓励参数为零。因此，非零参数对于回归问题来说才是重要的，这也正是我们称之为“变量选择”的原因。