概率预测与更多

在本书中，我们学习了生成预测的不同技术，包括一些经典方法、使用机器学习以及一些深度学习架构。但是，我们一直专注于一种典型的预测问题——为具有良好历史数据且无层次结构的连续时间序列生成点预测。我们这样做的原因是，这是你将面临的最常见的问题。但在本章中，我们将花时间讨论一些小众主题，这些主题虽然不那么流行，但同样重要。

在本章中，我们将重点讨论以下主题：

概率预测

到目前为止，我们一直在讨论将预测作为一个单一数字。我们一直在将深度学习模型投影到单一维度，或训练机器学习模型输出单一数字。随后，我们使用诸如均方损失之类的损失函数来训练模型。这种范式我们称之为点预测。但我们没有考虑一个重要方面。我们利用历史数据训练模型以做出最佳猜测。但模型对其预测有多确信？熟悉机器学习和分类问题的人会认识到，对于分类问题，除了得到样本属于哪个类别的预测外，我们还会得到模型不确定性的一些概念。但我们的预测是一个回归问题，我们无法免费获得不确定性。

但为什么量化不确定性在预测中很重要？任何预测都是为了某个目的而创建的，用于某些下游任务，这些任务会使用预测的信息。换句话说，利用我们生成的预测需要做出一些决策。而在做决策时，我们通常希望拥有尽可能多的信息。

让我们看一个例子来更清楚地说明这一点。你记录了最近5年的月度杂货消耗量，并利用本书中的技术创建了一个超级准确的预测和一个应用程序，告诉你每个月需要购买多少。你打开应用程序，它告诉你这个月需要购买两条面包。你前往超市，拿了两条面包回家。而在月底前一周，面包吃完了，剩下的时间你挨饿了。在饥饿的高峰期，你开始质疑自己的决策和预测。你分析数据找出哪里出错，意识到你每月的面包消耗量变化很大。有些月份你消耗了4条，而在其他月份只有1条。因此，这个预测很有可能在几个月里让你没有面包，而在另外几个月里有过多的面包。

然后，你读了这一章，将你的预测转换为概率预测，现在它告诉你，50%的情况下你下个月的面包消耗量是2条。但现在，应用程序增加了一个新功能，询问你更喜欢挨饿还是月底有剩余面包。因此，根据你对挨饿或省钱的倾向，你选择一个选项。假设你不想挨饿，但可以接受10%的情况下面包用完。一旦你将这个偏好输入应用程序，它会修正预测并告诉你应该购买3条面包，然后你再也不会挨饿了（同时因为你聪明了，从超市买了其他东西）。

利用预测中的不确定性并根据我们对风险的偏好进行修正是概率预测的主要效用之一。它还有助于让我们的预测对用户更加透明和可信。

现在，让我们快速看一下在使用学习模型的预测问题中可能遇到的不确定性类型。

预测不确定性类型

我们在第5章中看到，监督机器学习不过是学习一个函数 $\hat{y} = h(X, \phi)$，其中 h，连同 $\phi$，是我们学习的模型，而 X 是输入数据。因此，如果我们考虑不确定性的来源，它可以来自这两个组成部分中的任何一个。

我们所学习的模型 h 是对数据集 X 的一个近似，该数据集可能并未完全覆盖所有情况，因此系统中可能会引入一些不确定性。我们称此为认知不确定性。在机器学习背景下，当模型没有接触到足够的数据、模型本身不足以学习问题的复杂性，或者训练数据没有代表所有可能场景时，就可能出现认知不确定性。这也被称为系统性或可缩减的不确定性，因为这是总预测不确定性中可以通过更好的模型、更好的数据等主动减少的部分；换句话说，通过获得更多关于系统的知识。让我们看几个例子来使概念更清晰：

这类不确定性有一个好处，即我们可以通过收集更好的数据、训练更好的模型等方式主动减少它。

现在，总预测不确定性中还存在另一种不确定性——随机不确定性。这指的是数据中固有的、无法解释的随机性。这也被称为统计不确定性或不可约不确定性。尽管我们的宇宙看似是确定性的，但表面之下始终存在一层不确定性。

例如，天体运动可以精确计算（多亏了广义相对论和爱因斯坦），但一颗随机的小行星仍可能撞击任何天体，导致计算的轨迹发生改变。这种不可约且不可避免的不确定性被称为随机不确定性。来看几个例子：

现在了解了不同类型的不确定性以及为什么需要不确定性量化，让我们看看它在预测中的含义。

什么是概率预测和预测区间？

概率预测是指预测不仅给出单点预测，还捕捉了预测的不确定性。概率预测是一种预测未来事件或结果的方法，它提供一系列可能的值及其相关的概率或置信水平。这种方法捕捉了预测过程中的不确定性。

在计量经济学和经典时间序列领域，预测区间已经内置于公式中。这些方法的统计基础和强假设确保了模型的输出很容易以概率方式解释（只要你能满足这些模型规定的假设）。但在现代机器学习/深度学习领域，概率预测并非事后考虑。多个因素（如较少的刚性假设以及模型训练的方式）共同导致了这种困境。

有多种方法可以为预测添加概率维度，我们将在本章中介绍其中的几种。但在那之前，让我们先了解概率预测最有用的表现形式之一——预测区间。

预测区间是一个范围，未来观测值以指定概率落在此范围内。例如，如果[5,8]时间步的95%预测区间为5到8，那么95%的情况下实际值会落在5到8之间。我们以一个均值为$\mu_t$、方差为$\sigma_t^2$的正态分布为例，该分布作为时间步t的预测（我们将讨论的一种技术正好给出这个）。因此，在时间t，显著性水平（预测值落在区间外的概率）为$\alpha$的预测区间可以写为：

$$[\mu_t - z\sigma, \mu_t + z\sigma]$$

其中z是正态分布的z分数。

预测区间（PI）与置信区间（CI）

最容易混淆的概念之一是预测区间和置信区间。我们在此予以澄清。两者都是量化不确定性的方式，但在预测的语境中，它们服务于不同目的，且解释不同。置信区间也提供一个范围，但针对的是样本数据的总体参数（如均值），而预测区间侧重于为未来观测值提供范围。两者之间一个关键区别是，置信区间通常比预测区间窄，因为预测区间还考虑了新数据点的不确定性。而置信区间仅考虑模型参数的不确定性。因此，当需要给出未来观测值可能落入的范围（如下个月的销售额）时，我们使用预测区间。当需要为估计参数（如一年内的估计平均需求）提供范围时，我们使用置信区间。

在详细讨论之前，我们需要澄清几个术语和概念。

置信水平、错误率和分位数

在处理预测区间时，理解置信水平、错误率和分位数之间的关系至关重要。这些概念有助于确定未来观测值以一定概率落入的范围。

错误率（$\alpha$）是预测区间不包含未来观测值的允许概率。通常以百分比或0到1的小数表示。如果我们说$\alpha = 10\%$或$\alpha = 0.1$，意味着有10%的概率未来观测值不在预测区间内。

置信水平（$CL$）是错误率的补数，即预测区间包含未来观测值的概率。$CL = 1 - \alpha$。如果我们说错误率为10%，那么置信水平就是90%。

分位数是将数据分成等概率区间的点。简单来说，分位数表示低于该值的数据所占的百分比。例如，第5百分位数或0.05分位数表示5%的数据低于该点。因此，当我们无法基于分布假设解析地得到预测区间时，也可以用分位数来定义预测区间。

在图17.1中，我们展示了标准正态分布的预测区间。

图17.1：标准正态分布的预测区间

错误率、置信水平和分位数之间存在紧密联系。错误率和置信水平直接互补，可以互换使用，以定义我们希望预测区间具有的置信度，或者我们能够从预测区间接受的错误率。另一种理解方式是看曲线下的面积。在图17.1中，绿色阴影区域的面积表示置信水平，红色区域的面积表示错误率。

对于标准正态分布，我们可以直接通过解析公式得到预测区间：

$$[\mu - z_{\alpha/2} \cdot \sigma, \mu + z_{\alpha/2} \cdot \sigma]$$

其中 $\mu$ 是分布的均值，$\sigma$ 是分布的标准差，而 $z_{\alpha/2}$ 是标准正态分布中与期望置信水平对应的临界值，即 $1 - \alpha$。采用 $\alpha / 2$ 是因为我们允许错误率分布在两侧（即图17.1中曲线两侧的红色阴影区域）。

现在，让我们看看错误率和置信水平如何与分位数关联，因为如果我们不知道分布是什么（并且我们不想假设任何分布），我们就无法通过解析公式获得预测区间。在这种情况下，我们可以使用分位数来获得同样的结果。就像我们在解析公式中所做的那样，错误率 $\alpha$ 应平均分配到两侧。因此，预测区间为：

$$[q_{\alpha/2}, q_{1-(\alpha/2)}]$$

其中 $q_t$ 是第 tth 分位数。因此，根据分位数的定义，我们知道 $q_{\alpha/2}$ 以下有 $\alpha / 2$% 的数据，而 $q_{1-(\alpha/2)}$ 以上有 $\alpha / 2$% 的数据，从而使得区间外部的面积为 $\alpha$。

利用这一关系，我们可以从错误率转换为分位数，或者从置信水平转换为分位数。如果错误率为 $\alpha$，我们已经看到了相应的预测区间分位数。让我们再看一个从置信水平（以百分比表示）转换为分位数的快速公式：

$$[q_{50-(CL/2)}, q_{50+(CL/2)}]$$

在Python代码中，这很简单：

level = 95  # Confidence levels
qs = [50 - level / 2, 50 + level / 2] # Quantiles

现在，让我们看看如何衡量预测区间的质量。

测量预测区间的好坏

我们已经知道什么是概率预测和预测区间。但在探讨生成预测区间的技术之前，我们需要一种方法来衡量这类区间的好坏。诸如平均绝对误差或均方误差之类的标准指标不再适用，因为它们是点预测的衡量指标。

我们期望预测区间具备什么特性？如果我们有一个置信度为90%的预测区间，我们希望数据点至少有90%的时间落在该区间内。这很容易通过设置很宽的预测区间来实现，但这样一来预测区间就变得毫无用处。因此，我们希望预测区间尽可能窄，同时仍然满足90%的置信度标准。为了衡量这两个不同的方面，我们可以使用两个指标——覆盖率和预测区间的平均长度。

覆盖率是真值落在预测区间内的比例。数学上，如果我们用$[L_i, U_i]$表示每个观测点i的预测区间，用$y_i$表示真值，那么覆盖率可以定义为：

$$Coverage = \frac{1}{N} \sum_{i=1}^{N} \mathbb{I}(L_i \leq y_i \leq U_i)$$

其中$\mathbb{I}$是指示函数，若内部条件为真则等于1，否则为0，N是观测总数。覆盖率指标接近期望置信水平（例如，对于95%预测区间，接近95%）表明模型的置信度估计校准良好。

预测区间的平均长度是通过对所有观测点的预测区间长度取平均得到的。使用与上述相同的符号，可以数学表示为：

$$Average\ Length = \frac{1}{N} \sum_{i=1}^{N} (U_i - L_i)$$

该指标有助于理解区间覆盖率与其精确性之间的权衡。我们再来看看这两个指标的Python函数：

覆盖率：

import numpy as np
def coverage(y_true, lower_bounds, upper_bounds):
    """
    Calculate the coverage of prediction intervals.
    Parameters:
    y_true (array-like): True values.
    lower_bounds (array-like): Lower bounds of prediction intervals.
    upper_bounds (array-like): Upper bounds of prediction intervals.
    Returns:
    float: Coverage metric.
    """
    y_true = np.array(y_true)
    lower_bounds = np.array(lower_bounds)
    upper_bounds = np.array(upper_bounds)
    # Check if true values fall within the prediction intervals
    coverage = np.mean((y_true >= lower_bounds) & (y_true <= upper_bounds))
    return coverage

平均长度：

def average_length(lower_bounds, upper_bounds):
    """
    Calculate the average length of prediction intervals.
    Parameters:
    lower_bounds (array-like): Lower bounds of prediction intervals.
    upper_bounds (array-like): Upper bounds of prediction intervals.
    Returns:
    float: Average length of prediction intervals.
    """
    lower_bounds = np.array(lower_bounds)
    upper_bounds = np.array(upper_bounds)
    # Calculate the length of each prediction interval
    lengths = upper_bounds - lower_bounds
    # Calculate the average length
    average_length = np.mean(lengths)
    return average_length

这两个Python函数均可在src/utils/ts_utils.py中找到，我们将在本章中使用它们来衡量生成的预测区间的质量。

现在，我们来看可用于获取概率预测的不同技术，以及如何实际使用它们。

概率密度函数

这是概率预测中最常用的技术之一，特别是在深度学习领域，因为实现简单。时间t的预测值 $\hat{y}_i$ 可以看作概率分布 $p(\hat{y}_t)$ 的一个实现。我们不是直接估计 $\hat{y}_i$，而是估计 $p(\hat{y}_t)$。如果我们假设 $p(\hat{y}_t)$ 是参数化分布 $\mathcal{D}(\theta_t)$ 之一，参数为 $\theta_t$，那么我们可以直接估计参数 $\theta_t$，而不是 $\hat{y}_i$。

例如，如果我们假设预测值来自正态分布，我们可以建模为

$$\hat{y}_t \sim \mathcal{N}(\mu_t, \sigma_t)$$

因此，我们不是让模型输出 $\hat{y}_i$，而是让它输出 $\mu_t$ 和 $\sigma_t$。有了 $\mu_t$ 和 $\sigma_t$，我们可以轻松计算给定 $\alpha$ 的预测区间：

$$Lower\ Bound = \mu - Z \cdot \sigma$$

$$Upper\ Bound = \mu + Z \cdot \sigma$$

其中Z 是标准正态分布中对应于所需置信水平的临界值。对于90%置信水平，$Z \approx 1.645$。够简单吧？先别急！

现在我们要对分布的参数进行建模，那如何训练模型呢？我们仍然将实际点预测作为目标。在正态分布情况下，目标仍然是实际的 $y_t$，而不是均值和标准差。我们通过使用像对数似然这样的损失函数来解决这个问题，而不是均方误差之类的损失。

例如，假设我们有一组独立同分布的观测值（在我们的情况中就是目标），$y_1, y_2, \ldots, y_n$。利用假设分布的预测参数（$D(\theta)$），我们可以计算每个目标 $p(y_t)$ 的概率。独立同分布假设意味着每个样本相互独立。高中数学告诉我们，当两个独立事件发生时，我们可以将两个独立概率相乘来计算它们的联合概率。利用同样的逻辑，我们可以通过将所有单个概率相乘来计算所有n个独立同分布观测值的联合概率或似然（所有这些事件发生的概率）。

$$Likelihood = \prod_{i=1}^{n} p(y_i)$$

最大化似然有助于模型为每个样本学习正确的参数，使得在假设分布下的概率最大化。我们可以直观地这样理解：对于像正态分布这样的假设分布，最大化似然确保目标值落在由每个样本预测参数所定义的分布中心。

然而，这种操作在数值上并不稳定。由于概率 $0 \lt p \lt 1$ 相乘会使结果越来越小，很快导致数值下溢问题。因此，我们使用对数似然，它不过是似然的对数变换。这样做是因为：

$$Log\ Likelihood = \sum_{i=1}^{n} \log(p(y_t))$$

这种方法的主要缺点是它依赖于参数化的概率分布，其对数似然计算是可处理的。因此，我们被迫对输出做出假设，并提前选择一个可能合适的分布。

这是一把双刃剑。一方面，我们可以将一些领域知识注入问题中，并规范模型训练；但另一方面，如果我们不清楚选择正态分布是否正确，就可能导致模型受到不必要的约束。

许多流行的分布，如正态分布、泊松分布、负二项分布、指数分布、对数正态分布、Tweedie分布等，都可以用于通过这种技术生成概率预测。

现在，我们已经拥有了训练和学习模型的所有组件，并且可以使其预测完整的概率分布，而不是点预测。有了这些理论基础，让我们换个角度，看看如何使用这种技术。

基于概率密度函数（PDF）的预测——机器学习模型

我们在《第2部分：时间序列的机器学习》中已经看到如何使用标准机器学习模型进行预测。但所有这些都是点预测。我们能否轻松地使用PDF方法将其全部转换为概率预测？理论上可以，但实际并不容易。所有流行的机器学习模型实现，如 sci-kit learn、xgboost、lightgbm 等，都采用点预测范式。作为这些开源库的用户，我们很难调整和重写代码，使其将对数似然优化为回归损失。但别担心，并非不可能。NGBoost 是非常流行的梯度提升模型（如 xgboost 和 lightgbm）的近亲，其实现方式使其能够预测PDF，而不是点预测。

如果最终目标是获得预测区间，还有其他技术，如分位数预测或共形预测，对于本书中讨论的机器学习模型更为广泛适用（且推荐）。NGBoost的讨论是为了完整性，以及需要输出完整概率分布的情况。

我们不会深入探讨NGBoost是什么以及它与普通梯度提升模型有何不同，只需知道它是一个预测概率分布而非点预测的模型。

参考文献核实：

Duan等人提出的NGBoost研究论文在参考文献中的参考文献1中被引用。

延伸阅读提供了一个关于NGBoost的博客链接，该博客更深入地介绍了该模型是什么。

笔记本提示：

要跟随完整代码，请使用Chapter17文件夹中的名为01-NGBoost_prediction_intervals.ipynb的笔记本。

我们使用M4竞赛（包含10万个时间序列，参考文献5）中的8个时间序列样本进行概率预测。这些数据在线上很容易获取，下载脚本包含在笔记本中。我们使用这个比之前更简单的数据集，是为了避免用外生变量等复杂化叙述。以下是八个采样时间序列最后100个时间步的图。

图17.2：来自M4竞赛的8个采样时间序列的最后100个时间步。测试期用紫色虚线标出。

让我们使用mlforecast快速创建一些特征，以便将其转化为回归问题。我们在第6章有一个额外的笔记本展示了如何使用mlforecast作为本书仓库中包含的特征工程的替代方案。详细代码请参考完整笔记本，但现在假设我们有一个名为data的数据框，它包含了运行机器学习模型所需的所有特征。我们将其拆分为train和val，然后进一步拆分为X_train、y_train、X_val和y_val。现在，让我们看看如何训练一个假设输出为正态分布的模型。

from ngboost import NGBRegressor
from ngboost.distns import Normal
# Training the model
ngb = NGBRegressor(Dist=Normal).fit(X_train, Y_train)

NGBoost没有很多参数需要调整，因此不如其他梯度提升决策树（GBDT）灵活。而且它也不如其他GBDT快。这是一个仅在需要概率输出的特殊用例中使用的模型。

以下是NGBoost的一些参数：

现在我们已经训练了一个NGBoost模型，让我们看看如何使用它来生成预测和预测区间。

要获得点预测，语法与scikit-learn API完全相同。

y_pred = ngb.predict(X_val)

这只是一个包装方法，用于计算假定分布的位置参数。例如，对于正态分布，预测分布的平均值就是点预测。

现在要获取底层的概率预测以及随后的预测区间，我们需要使用一个不同的方法：

y_pred_dists = ngb.pred_dist(X_val)

y_pred_dists中的每个点都是一个完整的分布。如果我们想要查看前五个预测点，可以执行以下操作：

y_pred_dists[0:5].params

现在，为了获得预测区间，我们可以使用y_pred_dist并调用一个方法，给出我们期望的置信水平。这反过来会调用scipy分布（如scipy.stats.norm），它有一个方法interval来根据置信水平获取区间。

y_pred_lower, y_pred_upper = y_pred_dists.dist.interval(0.95)

现在，我们在y_pred周围有一个足够宽的窗口，以包含每个数据点预期的不确定性——预测区间。让我们看看预测结果和指标。

图17.3：来自NGBoost的预测及其预测区间

下面是针对这八个时间序列计算得到的指标（平均绝对误差、覆盖率及平均区间长度）：

图17.4：NGBoost的指标

虽然某些时间序列的预测区间看起来不错，但另一些（如时间序列H103）的预测区间似乎过窄，这一点同样在较低的覆盖率中体现。

基于PDF的预测——深度学习模型

与机器学习模型不同，将本书中学习过的所有深度学习模型转换为它们的PDF版本非常容易。还记得我们在开始实际应用之前的讨论吗？我们需要做的主要改动如下：

在深度学习范式中，这些改动非常简单，不是吗？在本书中我们学会使用的几乎所有深度学习模型中，最后都有一个线性投影层，将输出投影到所需维度。从概念上讲，只需改变线性投影即可使输出变为多个数值（假设分布的参数）。同样，改变损失函数也非常简单。值得注意的是，DeepAR （参考文献3）是一个著名的深度学习模型，它使用这种技术进行概率预测。

笔记本提示：

要跟随完整代码，请使用Chapter17文件夹中的名为02-NeuralForecast_prediction_intervals_PDF.ipynb的笔记本。

让我们看看如何在neuralforecast中实现这一点（我们在第16章中使用的库）。在我们的示例中，我们将采用一个简单的模型，如LSTM，但我们可以对任何模型做同样的操作，因为我们所做的只是将损失函数切换为DistributionLoss。

并且，在本例中，我们将使用M4竞赛数据集，该数据集可免费获取（下载数据集的代码包含在笔记本中）。

让我们从数据格式化满足neuralforecast在Y_train_df和Y_test_df中的要求开始。我们需要做的第一件事是导入必要的类。

from neuralforecast import NeuralForecast
from neuralforecast.models import LSTM
from neuralforecast.losses.pytorch import DistributionLoss

我们在本书中之前没有看到的唯一类是DistributionLoss。这是一个包装torch.distribution类并实现我们之前讨论的负对数似然损失的类。在撰写本书时，DistributionLoss类支持以下底层分布：

这些不同分布之间的选择完全取决于建模者，并且是模型中的一个关键假设。如果我们建模的输出预期服从正态分布，那么我们可以选择Normal。以下是DistributionLoss的主要参数：

现在，让我们设置一个预测期、所需的水平以及 LSTM 的几个超参数（我们选择了一些简单的小超参数以加快训练速度。在实际问题中，建议进行超参数搜索以找到最佳参数）。

horizon = 48
levels = [80, 90]
lstm_config = dict(input_size=3*horizon, encoder_hidden_size=8, decoder_hidden_size=8)

现在，我们需要定义将要使用的模型和 NeuralForecast 类。让我们定义两个模型——一个使用 Normal，另一个使用 StudentT。

models = [
    LSTM(
        h=horizon,
        loss=DistributionLoss(distribution="StudentT", level=levels),
        alias="LSTM_StudentT",
        **lstm_config
    ),
    LSTM(
        h=horizon,
        loss=DistributionLoss(distribution="Normal", level=levels),
        alias="LSTM_Normal",
        **lstm_config
    ),
]
# Setting freq=1 because the ds column is not date, but instead a sequentially increasing number
nf = NeuralForecast(models=models, freq=1)

注意，除了我们选择的分布损失函数外，语法与点预测完全相同。现在剩下的就是训练模型了。

nf.fit(df=Y_train_df)

模型训练完成后，我们可以使用 predict 方法进行预测。输出结果将包含我们定义的别名下的点预测，以及所有定义水平的高、低区间。

Y_hat_df = nf.predict()

现在我们有了概率预测，让我们看看它们并计算指标。

图17.5：LSTM 以 StudentT 分布作为输出的带预测区间的预测结果

图17.6：LSTM 以正态分布和 StudentT 分布输出的指标

如果我们将覆盖率与 NGBoost 的结果进行比较，可以看到深度学习方法提高了覆盖率，但在大多数情况下，区间也过宽（表现为较大的平均宽度）。

该方法最大的缺点是我们将输出限制为参数化分布之一。在许多实际案例中，数据可能不符合任何参数分布。

现在，让我们来看一种不需要假设任何参数分布，但仍能获取预测区间的方法。

分位数函数

如果我们的唯一目标是获得预测区间，我们也可以使用分位数来实现同样的目的。让我们从稍微不同的角度来审视PDF方法。在PDF方法中，我们在每个时间步输出一个完整的概率分布，并利用该分布得到分位数，即预测区间。尽管对于大多数参数分布，存在解析公式来计算分位数，但我们也可以数值求解。只需抽取N个样本（N足够大），然后计算这些样本的分位数。关键在于，即使我们拥有完整的概率分布，对于预测区间，我们需要的只是分位数。而给定N个样本计算分位数并不依赖于分布的类型。

那么，如果我们能训练模型直接预测指定的分位数，而不对潜在的概率分布做任何假设，会怎样呢？这正是分位数函数所做的。

在讨论分位数函数之前，我们先花一点时间了解累积分布函数（CDF）。这又是高中概率论的内容。简单来说，CDF返回某个随机变量X小于或等于某个值x的概率：

$$F(x) = P(X \leq x)$$

其中F是CDF。该函数接受输入x，并返回一个介于0和1之间的值。我们称这个值为$p$。

分位数函数是CDF的反函数。该函数告诉你，要使$F(x)$返回特定值$p$，x应取何值。

$$F^{-1}(p) = x$$

函数$F^{-1}$就是分位数函数。

从实现的角度看，对于能够进行多输出预测的模型（如深度学习模型），我们可以使用单一模型，通过改变输出层来预测所有我们想要的分位数。而对于限于单输出的模型（如机器学习模型），我们可以为每个分位数学习独立的分位数模型。

现在，和之前一样，我们不能使用像均方误差这样的点损失。在PDFs中，我们通过使用对数似然函数克服了这个问题。在这里，我们可以使用分位数损失或pinball损失。

参考文献核实：

提出分位数损失和回归的论文在参考文献的参考文献4中被引用。

分位数损失定义如下：

$$L_q(y_t, \hat{y}_t^q) = \begin{cases} (y_t - \hat{y}_t^q) \times q, & \text{if } y_t \geq \hat{y}_t^q \\ (\hat{y}_t^q - \hat{y}_t) \times (1 - q), & \text{if } y_t \lt \hat{y}_t^q \end{cases}$$

其中 $y_t$ 是目标值在时间 t 的目标值在时间 t，, $\hat{y}_t^q$ 是分位数预测，并且是分位数预测，q 是我们预测的分位数。这个公式看起来令人生畏，但请稍等片刻；它非常简单。让我们尝试获得关于损失的一些直觉。我们知道中位数是0.5分位数，这是一个集中趋势的度量。但如果我们要让预测接近75th百分位数或0.75分位数，那么我们必须促使模型高估，对吧？如果我们希望模型高估，我们需要在模型低估时给予更多惩罚。反之，如果我们要预测0.25分位数，我们需要低估。分位数损失正是这样做的：

不对称性来自项来自项 q 或 1 - 或 1 - q。另一项只是实际值和预测值之间的差异。

让我们尝试通过一个例子来理解这一点。假设我们有真实值 $y_t = 100$，并且我们想估计0.75分位数（

$$L_{q(100,110)} = (1 - 0.75) \times (110 - 100) = 2.5$$

$$L_{q(100,110)} = 0.75 \times (110 - 100) = 7.5$$

笔记本提示：

要跟随完整代码，请使用Chapter17文件夹中的名为03-Understanding_Quantile_Loss.ipynb的笔记本。

因此，通过改变q的值，我们可以使损失变得或多或少不对称，偏向任何一侧。下面的图17.7显示了q=0.5和q=0.75的损失曲线，并标出了这些示例预测。

图17.7：q=0.5和q=0.75的分位数损失曲线

我们可以看到，q=0.5的分位数损失是对称的，因为它代表中位数，而q=0.75的损失曲线是不对称的，对低估的惩罚远大于高估。

尽管公式具有分支结构，但在代码中实现时可以轻松避免取最大值操作。Python中的分位数损失如下所示：

def quantile_loss(q, y, y_hat_q):
    """ Calculate the quantile loss for a given quantile.
    Args:
    q (float): The quantile to be evaluated, e.g., 0.5 for median.
    y (float): The target value.
    y_hat_q (float): The quantile forecast.
    """
    error = y - y_hat_q
    return np.maximum(q * error, (q - 1) * error)

现在，让我们进入实际方面，学习如何将分位数损失与本书中涵盖的不同技术（机器学习和深度学习）结合使用。

使用分位数损失进行预测（机器学习）

常规机器学习模型通常只能建模一个输出。因此，我们需要为每个感兴趣的分位数训练一个不同的模型，使用分位数损失。也就是说，如果我们想预测某个问题的0.5、0.05和0.95分位数，我们必须训练三个独立的模型，每个分位数一个。

笔记本提示：

要跟随完整代码，请使用04-LightGBM_Prediction_Interval_Quantile_Loss.ipynb文件夹中的Chapter17笔记本。

让我们看看如何做到这一点。与PDF部分类似，我们使用mlforecast快速构建一个合成问题并创建一些特征。详细代码请参考完整笔记本，但现在，假设我们有一个数据框data，其中包含运行机器学习模型所需的所有特征。我们将其拆分为train和val，然后进一步拆分为X_train、y_train、X_val和y_val。

第一步是导入LGBMRegressor并设置一些参数和我们想要训练的分位数。

params = {
    'objective': 'quantile',
    'metric': 'quantile',
    'max_depth': 4,
    'num_leaves': 15,
    'learning_rate': 0.1,
    'n_estimators': 100,
    'boosting_type': 'gbdt'
}
# converting levels to quantiles
# For 90% Confidence - 0.05 for lower, 0.5 for median, and 0.95 for upper
quantiles = [0.5] + sum([level_to_quantiles(l) for l in levels], [])

这里需要注意的关键参数是objective设置为quantile，metric也设置为quantile。其余的LightGBM参数可以根据每个用例进行调整和微调。现在，让我们训练所有分位数模型。

# Training a model each for the quantiles
quantile_models = {}
for q in quantiles:
    model = LGBMRegressor(alpha=q, **params)
    model = model.fit(X_train, Y_train)
    quantile_models[q] = model

现在模型已经训练完成，我们可以从分位数模型中获得点预测和预测区间。

# Point Forecast using the 0.5 quantile model
y_pred = quantile_models[0.5].predict(X_val)
# Prediction Intervals using the 0.1 and 0.9 quantile models
y_pred_lower = quantile_models[0.1].predict(X_val)
y_pred_upper = quantile_models[0.9].predict(X_val)

现在，让我们看看预测结果和指标。

图17.8：使用LightGBM分位数回归的预测区间预测

图17.9：LightGBM分位数回归的指标

这种方法缺点是我们为每个分位数训练一个模型。这会很快变得难以管理。训练三个模型而不是一个模型会增加总训练时间。另一个问题是，由于这三个模型的训练方式不同，它们可能具有不同的属性，并且它们学习解决问题的方式也可能非常不同。由于这种不一致性，预测区间也可能出现一些问题。我们在图17.7的许多时间序列中可以清楚地看到锯齿状的预测区间；它们似乎与中位数预测脱节。这是深度学习世界中不存在的问题。

分位数损失预测（深度学习）

在深度学习模型中，我们使用共同的学习结构，并针对不同的分位数在共享投影上使用不同的线性投影。这确保了底层表示和学习在所有分位数上是共同的，并且可以产生更一致的分位数预测集。因此，对于本书中介绍的所有深度学习模型，我们可以通过做两件事将其转化为分位数预测模型：

  1. 我们不预测点预测（单个数字），而是预测概率分布的参数（一个或多个数字）。
  2. 不使用均方误差之类的点损失，而是使用对数似然之类的概率评分函数。

正如我们在PDF部分所做的那样，我们只需要切换neuralforecast中的损失函数。

笔记本提示：

要跟随完整代码，请使用Chapter17文件夹中的名为05-NeuralForecast_prediction_intervals_Quantile_Loss.ipynb的笔记本。

让我们看看如何在neuralforecast（我们在第16章中使用的库）中实现这一点。和之前一样，我们将采用一个简单的模型，如LSTM，以及M4竞赛数据集，但我们也可以对任何模型或任何数据集进行相同的操作，因为我们要做的只是将损失函数切换为MQLoss（多分位数损失）。

让我们从数据按照neuralforecast期望的格式化为Y_train_df和Y_test_df开始。我们需要做的第一件事是导入必要的类。

from neuralforecast import NeuralForecast
from neuralforecast.models import LSTM
from neuralforecast.losses.pytorch import MQLoss

我们之前没有见过的唯一类是MQLoss。这个类只是为我们刚刚讨论的多个分位数计算分位数损失（这是你通常训练模型的方式）。以下是MQLoss的主要参数：

现在，让我们设置一个时间范围、我们需要的水平以及LSTM的一些超参数。

horizon = 48
levels = [80, 90]
lstm_config = dict(input_size=3*horizon)

现在，我们需要定义将要使用的模型和NeuralForecast类。让我们只定义一个模型。

models = [LSTM(h=horizon, loss=MQLoss(level=levels), **lstm_config)]
# Setting freq=1 because the ds column is not date, but instead a sequentially increasing number
nf = NeuralForecast(models=models, freq=1)

请注意，除了我们选择的多分位数损失函数外，其语法与点预测完全相同。现在，剩下的就是训练模型了。

nf.fit(df=Y_train_df)

模型训练完成后，我们可以使用 predict 方法进行预测。该输出将在我们定义的别名下包含点预测，以及所有定义层级的高低区间。

Y_hat_df = nf.predict()

现在我们有了在不假设输出分布的情况下获取预测区间的方法，这在不确定潜在输出分布的实际场景中非常有价值。让我们看看生成的概率预测及其指标。

图17.10：使用分位数回归（深度学习）的预测区间预测

图17.11：分位数回归（深度学习）的指标

我们可以看到预测区间相互之间以及与中位数预测非常同步，不像使用LightGBM的独立模型那样脱节。这是因为所有分位数都在进行相同的学习，只是最终的投影头不同。总体而言，覆盖效果也更好。

还有另一种获取预测区间的方式，它极其简单，但由于PyTorch的典型使用方式而不易实现。接下来我们将介绍这种方法。

蒙特卡洛Dropout

Dropout是深度学习中一种非常受欢迎的正则化层。不深入细节，dropout正则化是指在训练过程中随机将网络的部分权重置为零（推理时关闭dropout）。直观上，这迫使模型不依赖少数权重，而是将权重的重要性分布到整个网络。从另一个角度看，我们是在应用一种正则化（非常类似于Ridge正则化），确保没有单个权重过高从而大幅度影响输出。

从技术上讲，除了将部分的随机权重置为零外，我们还会通过保留节点/权重的比例进行归一化，从而对每一层进行去偏。现在我们来形式化这一层。如果dropout概率为 $p$，并应用于中间激活，h，那么经过dropout后的激活 $h'$ 为：

$$h' = \begin{cases} 0 & \text{with } p \text{ probability} \\ \frac{h}{1 - p} & \text{otherwise} \end{cases}$$

为什么在应用 dropout 时需要缩放/归一化输出？直观的回答是为了确保在训练（dropout 激活时）和推理（dropout 关闭时）期间输出具有相同的尺度。详细的回答如下。

假设没有 dropout 时，节点的输出为 h。现在，使用 dropout 时，以概率 $p$ 该输出变为 0，以概率 $1 - p$ 该输出为 h。因此，该节点的期望值为：$\mathbb{E}[output] = (1 - p) \cdot h + p \cdot 0$。这意味着输出的平均值缩小了因子 $1 - p$，这是不希望的，因为这会改变训练和推理期间值的尺度。因此，解决方案是将在 dropout 中保留的节点的输出缩放 $1 - p$。

现在，我们了解了 dropout。但请记住，我们只在训练时将 dropout 用作正则化。2015 年，Yarin Gal 等人提出，传统的 dropout 还可以作为 高斯过程的贝叶斯近似。这个短语中包含了许多我们之前未接触过的术语。让我们短暂绕道，从高层次理解这些术语，并在 进一步阅读 中提供其他链接。

参考文献核实：

Yarin Gal 等人关于蒙特卡洛 dropout 的论文在 参考文献 中的参考文献 2 中被引用。

贝叶斯推断 是一种统计方法，它随着更多证据或信息的出现来更新假设的概率。它基于贝叶斯定理，该定理以数学形式表达了先验概率、似然和后验概率之间的关系。形式上，贝叶斯定理给出为：

$$p(\theta | D) = \frac{p(D | \theta) \cdot p(\theta)}{p(D)}$$

其中 $p(\theta | D)$ 是假设 $\theta$ 在给定数据 D 下的 后验概率，$p(D | \theta)$ 是 似然（在给定假设 $\theta$ 下观察到数据 D 的概率），$p(\theta)$ 是观察数据之前假设 $\theta$ 的 先验概率，而 $p(D)$ 是边际似然或 证据（在所有可能假设下观察到数据的总体概率）。

尽管它有一些特定的术语，但这非常直观，并提供了一种结构化的方法，在面对证据或数据时更新我们的先验信念。我们从先验分布 $p(\theta)$ 开始，它代表我们关于参数的初始信念。随着我们观察到数据 D，我们更新信念得到后验分布 $p(\theta | D)$。这个后验分布结合了先验信息和观测数据的似然，提供了关于参数的新的、更新的信念。进一步阅读 为感兴趣的人提供了更详细的解释。它还有来自 Seeing Theory 的一个页面，以直观的方式帮助可视化这些内容。

现在，让我们继续讨论 高斯过程（GP）。正如我们在 第5章 中所见，监督学习是学习一个函数 $\hat{y} = h(X, \phi)$，其中 $\hat{y}$ 是我们感兴趣预测的量，h 是我们学习的函数，X 是输入数据，而 $\phi$ 代表模型参数。因此，GP 将这个函数视为概率分布，并使用贝叶斯推断利用可用的训练数据更新函数的后验。这是一种截然不同的从数据中学习的方法，并且本质上是概率性的。

还剩一个术语，即近似。在许多情况下，贝叶斯模型中复杂的后验分布使其难以处理。因此，我们有一种称为 变分推断 的技术，其中我们使用一个已知的参数分布族 q，并找到最接近真实后验的成员。深入了解变分推断和 GP 的细节超出了本书的范围，但我为感兴趣的人在 进一步阅读 中添加了一些链接。

因此，回到 dropout，Yarin Gal 等人证明，在每个权重层之前定义了 dropout 层的神经网络实际上是一个 GP 的贝叶斯近似。所以，如果带有 dropout 的模型是 GP，而 GP 是函数上的后验，那么这应该给我们一个概率性输出，对吗？但在这里，我们没有像正态分布那样定义良好的参数化概率分布，以便解析地计算分布的属性（如均值、标准差或分位数）。我们该怎么做呢？

还记得分位数函数部分开头的讨论吗？我们说过，如果从某个分布中抽取了N个样本，并且这个N足够大，我们就可以近似该分布的性质。这有一个名称，叫做蒙特卡洛采样。蒙特卡洛采样是一种计算技术，通过从分布中生成大量随机样本来估计该分布的统计性质。将这个想法应用到启用丢弃法的神经网络中，我们可以使用蒙特卡洛采样来评估函数的后验概率分布的性质，这意味着在推理过程中需要保持丢弃法开启，并通过执行前向传播N次来从后验分布中采样。

理论上的合理性允许我们将丢弃法应用于任何神经网络，通过简单的操作获得不确定性估计。这种简单性难道不美吗？

所以，所有这些都归结为以下简单的步骤，以获得预测区间：

  1. 选择任何深度学习架构。
  2. 在每个主要操作之前插入丢弃层，并将其设置为$p$，其中$p \gt 0$。
  3. 训练模型后，在启用丢弃法的情况下，执行N次前向传播。
  4. 利用这N个样本，估计中位数（用于点预测），以及对应于定义置信水平的预测区间的分位数。

这听起来足够简单，但仅因一个原因而变得复杂。PyTorch和Tensorflow的设计使得在推理过程中丢弃层被关闭。在PyTorch中，我们可以通过执行model.train()或model.eval()来指示模型处于训练阶段还是推理阶段。而大多数流行的封装PyTorch以简化并自动化训练的实现（如PyTorch Lightning）会在预测前在后台执行这个model.eval()步骤。因此，当使用像neuralforecast这样的库时（它在后台使用PyTorch Lightning），在预测期间开启丢弃法并不容易。

在我们学习如何为neuralforecast模型实现MC Dropout之前，我们先绕个小弯，学习如何在neuralforecast中定义自定义模型。当你想要针对自己的用例调整任何模型时，这很有用。我们在这里这样做有两个原因：

  1. 我想向你展示如何在neuralforecast中定义一个新模型。
  2. 我想要一个适合MC Dropout技术的模型，即模型需要在每个权重/层之前都有丢弃层。

在neuralforecast中创建自定义模型

现在是时候找点乐子，定义一个能与neuralforecast配合使用的自定义PyTorch模型。为了简洁起见，我们使用一个对D-Linear进行小修改的模型，该模型在第16章中学习过。除了线性趋势和季节性外，我们还增加了一个非线性趋势组件。我们也给它起个古怪的名字——D-NonLinear。其架构如下：

图17.12：D-NonLinear 模型架构

现在，我们来了解如何编写一个能与 neuralforecast 配合使用的模型。neuralforecast 中的所有模型都继承自三个基类之一——BaseWindows、BaseRecurrent 或 BaseMultivariate。文档清楚地说明了 BaseWindows 的用途，这正是我们用例所需要的。我们需要在训练时从时间序列中采样窗口。

我们还需要记住一点：neuralforecast 在底层使用 PyTorch Lightning 进行训练。这个链接有关于如何为 neuralforecast 定义新模型的更多详细信息：https://nixtlaverse.nixtla.io/neuralforecast/docs/tutorials/adding_models.html。

如果你不熟悉 面向对象编程（OOP）和 继承，那么你可能很难理解我们在这里做什么。继承允许子类继承父类中定义的所有属性和方法。这使得开发人员可以在基类中定义通用功能，然后继承该类以获得所有功能，并在此基础上添加任何特定功能。强烈建议你理解继承，不仅是为了这个示例，也是为了成为更好的开发者。互联网上有数百个教程，我在这里引用一个：https://ioflood.com/blog/python-inheritance。

模型的完整代码可以在 src/dl/nf_models.py 中找到，但我们在这里将查看模型定义的关键部分。

让我们从定义 __init__ 函数开始（这里只包含相关部分；完整类定义请参考 Python 文件）。

class DNonLinear(BaseWindows):
    def __init__(
        self,
        # Inherited hyperparameters with no defaults
        h,
        input_size,
        # Model specific hyperparameters
        # Window over which the moving average operates for trend extraction
        moving_avg_window=3,
        dropout=0.1,
        # Inhereted hyperparameters with defaults
        ...
        **trainer_kwargs,
    ):
        super(DropoutDNonLinear, self).__init__(
            h=h,
            ...
            **trainer_kwargs,
        )
        # Model specific hyperparameters
        self.moving_avg_window = moving_avg_window
        self.dropout = dropout
        # Model initialization to follow

我们现在已经定义了 __init__ 函数的一部分。现在，让我们初始化该方法其余部分所需的不同层。我们有一个序列分解层，它使用移动平均将输入拆分为趋势和季节性分量；一个线性趋势预测器和季节性预测器，它们接收线性趋势和季节性并将其投影到未来；以及一个非线性预测器，它接收原始输入并将其投影到未来。

# Defining a decomposition Layer
        self.decomp = SeriesDecomp(self.moving_avg_window)
        # Defining a non-linear trend predictor with dropout
        self.non_linear_block = nn.Sequential(
            nn.Dropout(self.dropout),
            nn.Linear(self.input_size, 100),
            nn.ReLU(),
            nn.Dropout(self.dropout),
            nn.Linear(100, 100),
            nn.ReLU(),
            nn.Dropout(self.dropout),
            nn.Linear(100, self.h),
        )
        # Defining a linear trend predictor with dropout
        self.linear_trend = nn.Sequential(
            nn.Dropout(self.dropout),
            nn.Linear(self.input_size, self.h),
        )
        # Defining a seasonality predictor with dropout
        self.seasonality = nn.Sequential(
            nn.Dropout(self.dropout),
            nn.Linear(self.input_size, self.h),
        )

现在，让我们定义 forward 方法。forward 方法应该只有一个参数，即一个包含不同输入的字典：

对于此用例，我们只需要使用用例，我们只需要 insample_y，因为我们的模型不使用任何其他信息。因此，这就是forward 方法实现：方法实现：

def forward(self, windows_batch):
        # Parse windows_batch
        insample_y = windows_batch[
            "insample_y"
        ].clone()  # --> (batch_size, input_size)
        seasonal_init, trend_init = self.decomp(
            insample_y
        )  # --> (batch_size, input_size)
        # Non-linear block
        non_linear_part = self.non_linear_block(
            insample_y
        )  # --> (batch_size, horizon)
        # Linear trend block
        trend_part = self.linear_trend(trend_init)  # --> (batch_size, horizon)
        # Seasonality block
        seasonal_part = self.seasonality(
            seasonal_init
        )  # --> (batch_size, horizon)
        # Combine the components
        forecast = (
            trend_part + seasonal_part + non_linear_part
        )  # --> (batch_size, horizon)
        # Map the forecast to the domain of the target
        forecast = self.loss.domain_map(forecast)
        return forecast

代码非常直接。我们唯一需要确保与 neuralforecast 模型对齐的事情是从输入字典中获取我们需要的数据，并在最后调用self.loss.domain_map，以便根据损失将其映射到正确的输出大小。现在，该模型将像 neuralforecast 库中的任何其他模型一样工作。

现在，让我们回到 MC Dropout 及其实现。

使用 MC Dropout 进行预测（neuralforecast）

我们知道，在 neuralforecast 和 PyTorch Lightning 这样的框架中实现 MC Dropout 并不容易，但困难并不意味着应该放弃。关键是在预测期间启用 dropout，并进行多次采样。如果你编写自己的 PyTorch 训练代码，只需在预测前不要调用 model.eval() 即可。不过，更佳的做法是只将 dropout 层设为训练模式，而不是让整个模型处于训练模式，因为 batch normalization 等层在推理期间的行为也会随之改变。下面介绍一种方便地将所有 dropout 层设为训练模式的方法。

def enable_dropout(model):
    """Function to enable the dropout layers during test-time"""
    for m in model.modules():
        if m.__class__.__name__.startswith("Dropout"):
            m.train()

对于 neuralforecast，我们整理了一个配方，可以使用 MC Dropout 对他们的任何模型（只要它们有足够的 dropout）。现在，我们将使用我们刚刚定义的自定义DNonLinear 模型。请注意，定义中的每个组件都以 dropout 层开头，这样我们就可以毫无顾虑地应用 MC Dropout。

笔记本提示：

要跟随完整代码，请使用Chapter17文件夹中的名为06-Prediction_Intervals_MCDropout.ipynb的笔记本。

如果你还记得，我们之前在第13章中使用了 PyTorch Lightning，并解释过它本质上是以特定形式组织的标准 PyTorch 代码，例如 training_step、validation_step、predict_step 和 configure_optimizers。如需复习，请返回第13章和进一步阅读部分，了解如何从标准 PyTorch 迁移到 PyTorch Lightning。由于 neuralforecast 的后端已经使用 PyTorch Lightning，我们继承的 BaseWindows 本身就是一个 PyTorch Lightning 模型。这一点很重要，因为我们需要修改其 predict_step 方法来实现 MC Dropout。

使用我们过去用于继承的同一继承，我们过去用于继承 BaseWindows，我们可以继承，我们可以继承 DNonLinear 类，我们之前定义的类，并做一些修改，使其成为 MC Dropout 模型。为此，我们需要重新定义的是predict_step 方法。该方法。 predict_step 方法是方法是 PyTorch Lightning 每次需要获取批次预测时调用的方法。因此，我们不应直接取预测结果，而需要保持 dropout 启用，从N 次前向传播中采样N 个样本，计算预测区间和中位数（点预测），然后返回。

class MCDropoutDNonLinear(DNonLinear):
    def predict_step(self, batch, batch_idx):
        enable_dropout(self)
        pred_samples = []
        # num_samples and levels will be saved to the model in MCNeuralForecast predict method
        for i in range(self.num_samples):
            y_hat = super().predict_step(batch, batch_idx)
            pred_samples.append(y_hat)
        # Stack the samples
        pred_samples = torch.stack(pred_samples, dim=0)
        # Calculate the median and the quantiles
        y_hat = [pred_samples.quantile(0.5, dim=0)]
        if self.levels is not None:
            for l in self.levels:
                lo, hi = level_to_quantiles(l)
                y_hat_lo = pred_samples.quantile(lo, dim=0)
                y_hat_hi = pred_samples.quantile(hi, dim=0)
                y_hat.extend([y_hat_lo, y_hat_hi])
        # Stack the results
        y_hat = torch.stack(y_hat, dim=-1)
        return y_hat

我们完成了吗？还没有。还有一件小事要做。在第16章中，我们使用了一个名为NeuralForecast 的类来拟合和预测neuralforecast 模型。这是一个包装器类，负责在调用底层模型之前以正确的方式准备输入和输出的繁重工作。这个类必须知道我们已经修改了predict_step，因此我们需要在那里做一个小的改动。该解决方案更像是一个 hack 而不是一个有原则的编辑方式，但如果 hack 能达到目的，它也一样好。我一直在审视实现，以找出最好的方式来 hackNeuralForecast，以启用我们的 MCDropout 推断。没有简短的方式解释这个 hack，只需理解我滥用了neuralforecast 根据不同的损失灵活产生点预测和预测区间的方式。所以，这里是重新定义的NeuralForecast 类，其中包含一个 hack 在predict 方法中。

class MCNeuralForecast(NeuralForecast):
    def __init__(self, num_samples, levels=None, **kwargs):
        super().__init__(**kwargs)
        self.num_samples = num_samples
        self.levels = levels
    def predict(
        self,
        df=None,
        static_df=None,
        futr_df=None,
        sort_df=True,
        verbose=False,
        engine=None,
        **data_kwargs,
    ):
        # Adding model columns to loss output names
        # Necessary hack to get the quantiles and format it correctly
        for model in self.models:
            model.loss.output_names = ["-median"]
            for l in list(self.levels):
                model.loss.output_names.append(f"-lo-{l}")
                model.loss.output_names.append(f"-hi-{l}")
            # Setting the number of samples and levels in the model
            model.num_samples = self.num_samples
            model.levels = self.levels
        return super().predict(
            df, static_df, futr_df, sort_df, verbose, engine, **data_kwargs
        )

就这样。我们已经成功“黑掉”了库，让它为我们所用。这不仅是关于MC Dropout的教程，也是关于如何黑掉库以实现你期望功能的教程。值得注意的是，这并不会让你成为“黑客”，所以在你更新领英头衔之前，请先三思。

现在，开始训练模型。这与在 neuralforecast 中训练其他模型基本相同，但你需要使用我们新定义的 MCNeuralForecast 类，而不是 NeuralForecast 类。

horizon = len(Y_test_df.ds.unique()) # 48
levels = [80, 90]
model = MCDropoutDNonLinear (
    h= horizon,
    input_size=WINDOW,
    moving_avg_window=horizon*3,
    dropout=0.1,
    max_steps=500,
    early_stop_patience_steps=5,
)
mcnf = MCNeuralForecast(models=[model], freq=1, num_samples=100, levels=levels)
mcnf.fit(Y_train_df, val_size= horizon, verbose=True)

训练完成后，我们可以这样生成预测：

Y_hat_df = mcnf.predict()
Y_hat_df = Y_hat_df.reset_index()

输出与 neuralforecast 中的其他模型完全一致，预测区间格式化为 &lt;ModelName&gt;-lo-&lt;level&gt; 和 &lt;ModelName&gt;-hi-&lt;level&gt;。点预测可以在 &lt;ModelName&gt;-median 中找到。在这种情况下，&lt;ModelName&gt; 将是 MCDropoutDNonLinear。

让我们看看预测图和指标。

图17.13：使用MC Dropout的预测区间

图17.14：MC Dropout的指标

我们的预测模型在数据上表现尚可；不算出类拔萃，但还可以。如果进行消融研究，我们甚至可能发现添加的非线性组件毫无作用。但只要我们在过程中感到乐趣并学到东西，我就很开心。现在，看看预测区间。它们并不平滑，而且看起来有不少“噪声”，对吧？这是因为方法本身固有的随机性，也可能是学习不充分造成的。当我们进行MC Dropout时，本质上依赖于 N 个子模型或子网络，并基于这些 N 个预测计算分位数。可能其中一些子网络学习得不够好，它们的输出会扭曲分位数，从而影响预测区间。

MC Dropout方法受到许多批评。贝叶斯学派的许多人并不认为MC Dropout是贝叶斯方法，并且认为所提出的变分近似是一个非常粗糙的近似，以至于我们不能将其衡量的东西称为贝叶斯不确定性。有一篇Loic Le Folgoc等人未发表的Arxiv论文，题为“Is MC Dropout Bayesian?”（参考文献6），声称它不是。但这并不能改变MC Dropout是一种廉价的不确定性量化方法的事实。然而，在医学研究等不确定性量化至关重要的领域使用时，我们可能需要采取更严谨的方法。

我们还可以注意到，所有时间序列的覆盖率都很差。这也是不同研究中普遍存在的现象。2023年，Nicolas Dewolf等人发表了一项研究，比较了回归问题中不同的不确定性量化方法（参考文献7）。他们发现MC Dropout在覆盖率和平均长度方面的表现是最差的之一，这进一步证实了MC Dropout是一个非常粗糙的不确定性近似。

现在，让我们看看另一种概率预测技术，它承诺具有完美覆盖率的理论保证，并且在过去几年中变得非常流行。

共形预测

如果我告诉你有一种技术可以生成统计上保证完美覆盖的预测区间，适用于任何模型，并且不需要我们对输出分布做出任何假设，你会怎么想？共形预测正是如此。共形预测是一种通过估计模型的不确定性来帮助机器学习模型做出可靠预测的方法。它为任何机器学习模型提供稳健、统计上有效的不确定性度量，确保在关键应用中做出可靠且值得信赖的预测。

尽管它早在2005年就由Vladmir Vovk提出(参考文献 8)，但在过去几年才引起关注。让我们先通过一个分类示例理解共形预测的基本原理，然后看看如何将其应用于回归和时间序列示例。

共形预测用于分类

让我们从一个训练好的模型开始，$\hat{f}$，它输出估计的概率（softmax得分）用于K个输出类别（$\hat{f}(x) \in [0,1]^K$）。这个模型是什么并不重要；它可以是机器学习模型、深度学习模型，甚至是基于规则的模型。我们有训练数据$(X_{train}, Y_{train})$和测试数据$(X_{test}, Y_{test})$。现在，我们需要少量额外的数据（除了训练和测试数据之外），称为校准数据，$(X_{calib}, Y_{calib}) = \{(x_1, y_1), ..., (x_n, y_n)\}$。那么，我们想要从中得到什么？利用$\hat{f}$和$X_{calib}$，我们希望创建一个可能标签的预测集$\mathcal{C}(X_{test})$，确保测试数据点属于该集的概率几乎正好等于用户定义的错误率$\alpha$（在整个讨论中我们将讨论错误率。10%的错误率对应90%的置信水平）。这正是共形预测所保证的。这被称为边际覆盖保证，可以更正式地写成：

$$1 - \alpha \leq \mathbb{P}(Y_{test} \in \mathcal{C}(X_{test})) \leq 1 - \alpha + \frac{1}{n+1}$$

自然，你可能会想到这个问题。这个$\frac{1}{n+1}$是什么？这个项表示覆盖保证来自大小为n的有限样本。由于n在分母中，我们知道随着n增加，这项会越来越小。扩展到极限情况，如果$n = \infty$，这一项将为零，覆盖将恰好是$1 - \alpha$。

我们知道想要什么，但如何实现呢？共形预测的核心思想非常简单，可以分为四个步骤：

  1. 利用训练好的模型识别一种启发式的不确定性度量。在我们的分类示例中，这可以是softmax得分。
  2. 定义一个得分函数$s(\hat{y}, y) \in \mathbb{R}$，也称为非一致性分数。这是一个接收预测值$\hat{y}$和实际值$y$并给出编码它们之间不一致程度的分数。分数越高，不一致越大。在分类示例中，这可以是像$1 - \hat{f}(x)_{Y_i}$这样简单的东西。通俗地说，这表示取正确类别的softmax得分，然后计算$1 - score$。
  3. 计算$\hat{q}$作为校准分数的$\frac{(1 - \alpha)(n + 1)}{n}$分位数。我们使用校准数据和得分函数计算校准分数，然后计算这些数据的分位数。$\frac{n + 1}{n}$中的分位数计算同样源自有限样本修正。当$n$趋于无穷时，该项趋于零。
  4. 使用这个分位数为新样本形成预测集：$\mathcal{C}(x) = \{y : s(x, y) \leq \hat{q}\}$。这意味着从输出集中选择所有得分（根据得分函数）大于阈值$\hat{q}$的项。

这种简单技术将为我们提供保证满足边际覆盖率的预测集，无论使用什么模型或数据分布如何。让我们看看使用Python代码并假设我们讨论的模型是一个scikit-learn分类器，这有多简单。

我们有一个训练好的模型model、校准数据X_calib和测试数据X_test。完整代码和一些可视化内容，请查看笔记本。

笔记本提示：

要跟随完整代码，请使用Chapter17文件夹中的名为07-Understanding_Conformal_Prediction.ipynb的笔记本。

# 1: Get conformal scores
n = calib_y.shape[0]
cal_smx = model.predict_proba(calib_x) # shape (n, n_classes)
# scores from the softmax score for the correct class
cal_scores = 1 - cal_smx[np.arange(n), calib_y] # shape (n,)
# 2: Get adjusted quantile
alpha = 0.1  # Confidence level (1 - alpha)
q_level = np.ceil((n + 1) * (1 - alpha)) / n
qhat = np.quantile(cal_scores, q_level, method='higher')
# 3: Form prediction sets
val_smx = model.predict_proba(test_x)
prediction_sets = val_smx >= (1 - qhat)

现在，让我们考虑预测集，$\mathcal{C}$。我们一直将其定义为针对分类场景的离散类别的集合值。该集合根据初始启发式不确定性估计的置信度而变大或变小。

回归中的保形预测

让我们将预测集的概念扩展到回归。在回归中，输出空间是连续的而非离散的，我们的目标是构建连续的预测集，通常是$\mathbb{R}$中的一个连续区间。其思想是保持相同的覆盖原则：预测区间应以高概率包含真实值。因此，之前看到的预测集$\mathcal{C}$在回归上下文中也是预测区间。但随着预测集解释的变化，我们还需要改变计算非一致性得分的得分函数。一个常用的得分函数是到条件均值的距离，$s(x, y) = |y - \mu(x)|$（参考文献10）。当我们有一个训练好的模型$\hat{f}$时，我们可以将模型输出视为条件均值，这将使每个点的绝对残差成为该得分。

$$s(x, y) = |y - \hat{f}(x)|$$

注意，该得分满足条件。偏差越大，“启发式”不确定性度量越大。其余步骤几乎相同——计算分位数$\hat{q}$，并形成预测区间$\mathcal{C}(x) = [\hat{f}(x) - \hat{q}, \hat{f}(x) + \hat{q}]$。

让我们看看Python代码如何变化（完整代码在笔记本中）。

# 1: Get conformal scores
calib_preds = model.predict(calib_x)
cal_scores = np.abs(calib_y - calib_preds)
# 2: Get adjusted quantile
qhat = … # Exactly the same as classification
# 3: Form prediction intervals
test_preds = model.predict(test_x)
lower_bounds = test_preds - qhat
upper_bounds = test_preds + qhat

就这么简单。我们可以检查覆盖率，会看到它大于90%，这是我们用$\alpha = 0.1$定义的错误率。

从业者笔记

在许多使用场景中，我们将针对多个实体或我们关心的组训练单一模型。例如，对于我们在第10章中讨论的全局预测模型，我们为多个时间序列使用单一回归模型。在这种情况下，我们也可以分别对每个时间序列或时间序列组运行共形预测，以更适应于该子集的误差。这将允许在组或时间序列层面上提供覆盖保证。

但你可能已经注意到了一些问题。在这种方法中，区间宽度在整个过程中是相同的。但我们期望当模型更自信时区间更窄，不自信时区间更宽（从现在起我们称之为自适应预测区间）。让我们看看另一种具有此特性的技术。

共形分位数回归

我们了解了分位数回归作为概率预测（或一般情况下的回归）的方法之一。分位数回归的强大之处在于它不需要我们对底层输出分布有任何先验假设。但它不具备共形预测所提供的覆盖保证。2019年，Yaniv Roano等人（参考文献11）将两者的优势结合到了共形化分位数回归（CQR）中。他们提出了一种方法，即取分位数预测并进行共形化，使其具有共形预测所保证的覆盖保证。

在这种情况下，所使用的模型有一个限制。它应该是一个输出分位数预测的模型。而且，从我们之前的讨论中，我们知道如果误差率为$\alpha$，那么我们预测区间所需的分位数是$\hat{y}_t^{\alpha/2}$和$\hat{y}_t^{1 - (\alpha/2)}$。

因此，分位数模型预测$\hat{y}_t^{\alpha/2}$和$\hat{y}_t^{1 - (\alpha/2)}$。根据定义，如果$\hat{y}_t^{\alpha/2}$和$\hat{y}_t^{1 - (\alpha/2)}$是真实分位数的正确估计，那么仅使用分位数回归就会具有完美的覆盖。但模型拟合可能并不完美，这将导致覆盖不佳。我们将使用共形预测基于校准数据来校正分位数，从而获得共形预测所承诺的完美覆盖。

Yaniv Roano等人提出了一种用于分位数回归的新的非一致性得分函数。

$$s(x, y) = \max \{\hat{y}_t^{\alpha/2} - y, y - \hat{y}_t^{1 - (\alpha/2)}\}$$

让我们稍作停顿，使用图17.15中的示意图来探索这个得分函数。

图17.15：共形化分位数回归的得分函数说明

max运算符中有两项。如果真实值$y$位于两个分位数$\hat{y}_t^{\alpha/2}$和$\hat{y}_t^{1 - (\alpha/2)}$之间，则内部的两项都为负，并且将是到最近预测区间的距离（参见图17.15中的点B和C）。

现在，让我们看看点A（3），它位于较高分位数之上。分位数为[1, 2]。max运算符中的两项将是：

$$\hat{y}_t^{\alpha/2} - y = 1 - 3 = -2$$

$$y - \hat{y}_t^{1-(\alpha/2)} = 3 - 2 = 1$$

由于max算子，评分函数值为3。现在，我们看点D（1.6），它低于下分位数（分位数为[3.5, 2]）。这使得评分函数值为max{0.4, -1.9}，即0.4。因此，max算子确保如果值落在分位数之外，评分是正的。

因此，我们所拥有的是一个评分函数，它为实际值落在区间之外的点分配正值，为落在区间内的点分配负值。对于落在区间外的点，评分函数的构造方式会选择更严重的误差。这满足我们对评分函数的要求。更大的评分表示更大的不确定性，同时也编码了一种启发式的不确定性概念。

现在我们有了评分，剩下的步骤几乎相同：

  1. 计算 $\hat{q}$ 作为校准评分的 $\frac{(1 - \alpha)(n + 1)}{n}$ 分位数。
  2. 使用这个分位数来为新样本构建预测集：$\mathcal{C}(X_{test}) = [\hat{y}_t^{\alpha/2} - \hat{q}, \hat{y}_t^{1-(\alpha/2)} + \hat{q}]$，即我们将现有的分位数扩大 $\hat{q}$，并得到覆盖保证。

分位数回归是获得自适应预测区间的较好方法之一。现在我们来看这样一种技术。

共形不确定性估计

如果我们深入思考一下上一节中的共形分位数回归（Conformalized Quantile Regression），就会发现底层分位数回归所做的是捕捉预测中每个点的不确定性估计。然后我们将这些估计进行共形化以获得更好的覆盖。如果我们能捕捉到这种不确定性估计，即 $u(x)$，我们就有了将其共形化以获得更好的自适应预测区间的希望。

假设我们有一个训练好的模型 $\hat{f}$ 和一个不确定性标量 $u(x)$，当不确定性高时其值高，反之亦然。我们可以将非一致性评分定义为：

$$s(x, y) = \frac{|y - \hat{f}(x)|}{u(x)}$$

该分数的自然解释是，我们正在对标准 $|y - \hat{f}(x)|$ 乘以一个修正因子。一旦我们有了这个新分数，后续过程完全相同——从分数中取 $\hat{q}$，并形成预测区间 $\mathcal{C}(x) = [\hat{f}(x) - u(x) \cdot \hat{q}, \hat{f}(x) + u(x) \cdot \hat{q}]$。

那么，有哪些方法可以捕捉这种不确定性？（注意，这种不确定性度量应在数据点层面进行捕捉。）

尽管上述列表并不详尽，但它表明我们可以将共形预测应用于几乎任何不确定性估计（包括本章已经见过的那些）。这使得共形预测范式成为一个非常灵活的工具包，能够为各种问题提供覆盖保证。然而，即使具有这种灵活性，仍存在一些情况会破坏该框架所承诺的覆盖保证。对于我们的时间序列场景，理解这一点非常重要。

共形预测中的可交换性与时间序列预测

可交换性是共形预测中的一个基本假设。可交换性意味着数据点是独立同分布的，并且当数据点的顺序改变时，它们的联合概率分布不变。这一概念保证了过去的数据点能够可靠地预测未来的数据点。

想象一家巧克力工厂生产重量一致的巧克力。如果生产过程高度受控，且条件和原料相同，那么巧克力的重量是可交换的，因为生产顺序不影响重量。你可以抽取100块巧克力，测量其重量，并基于与预测重量的偏差计算非一致性分数。利用这些分数，你可以为未来的巧克力形成预测区间。由于巧克力是可交换的，样本分布代表了未来分布，因此预测区间是可靠的。

然而，如果生产过程随时间变化——例如由于机械磨损或不同批次的原料——那么重量就变得不可交换。生产顺序影响重量，使得样本分布无法代表未来的重量，从而导致不可靠的预测区间。

在时间序列数据中，观测值通常依赖于之前的观测值，这违反了可交换性假设。例如，在销售预测中，今天的销售可能因趋势影响明天的销售，或者一年前的销售可能因季节性效应影响明天的销售。这种依赖性意味着过去数据的分布并不能准确代表未来数据的分布。

但这对我们意味着什么？最明显的答案是，我们的覆盖保证会受到影响。但我们仍然可以将这些技术应用于时间序列吗？当然可以。经验上，学界发现这个框架同样适用于时间序列数据，只是覆盖保证会有所损失。在大多数实际应用中，对时间序列数据使用常规共形预测应该没有问题。2023年，Barber等人（参考文献12）研究了这个问题，并为非可交换数据（如时间序列）推导了理论覆盖保证。他们将覆盖缺口定义为期望覆盖率（$1 - \alpha$）与实际覆盖率之差，并推导了该缺口的上界，以显示分数可交换性假设被违反的程度。对于这个界限，他们考虑了使用原始模型的校准数据分数$s(z)$，以及另一个模型，该模型在相同数据上训练，但将训练数据中的一个随机选择的数据点与测试数据点$s(z^i)$交换。

该界限被证明与$d_{TV}(s(z), s(z^i))$成正比，即这两个分数之间的分布距离。在我们使用的大多数算法中，交换一个数据点不会显著改变模型，因此我们仍然可以使用共形预测处理时间序列数据，且覆盖损失最小。

但另一方面，如果我们希望预测区间非常准确，或者我们使用的模型特别容易受到数据点交换的影响，那么就需要一些技术来克服由于分布变化导致的性能下降。处理这个问题的方法有很多，在撰写本书时，这是一个活跃的研究领域。这里有两个非常简单的方法值得一提。

加权共形预测

假设一个时间序列的数据分布存在缓慢变化，$(x_1, x_2, ..., x_N)$，我们使用校准集$X_c = (x_{N-k}, ..., x_N)$，取最后k个时间步。并且我们感兴趣的是预测测试集$X_{test} = (x_{N+1}, ..., x_{N+H})$，其中H是预测时域。

因此，有理由认为$X_c$中最新的时间步最接近测试时间段内观察到的值分布。那么，如果我们为校准数据中的非一致性分数分配权重，使得最近的时间步获得更高的权重，并计算加权分位数而不是常规分位数，会怎样？显然，这是一个非常好的主意，并且有坚实的理论支持。

使用最近时间来加权校准数据只是我们利用权重处理分布偏移的一种方式。更一般地，任何权重调度$w_1, ...w_k, w_i \in [0, 1]$都可以在这里使用。也许对于强季节性的时间序列，使用季节性周期来定义权重是有意义的，或者可能存在其他已知的标准，使得校准数据的不同实例对未来预测的相关性不同。

在进入具体机制之前，我们需要了解什么是加权分位数。如果你已经熟悉这个概念，可以跳过。如果你需要一些直觉理解，我强烈建议你查看本章文件夹中名为08-Quantiles_and_Weighted_Quantiles.ipynb的笔记本。

现在，让我们回到加权共形预测方法。

因此，正如我们之前讨论的，对于任何归一化的权重调度$w_i \in [0,1] \text{ and } \sum_{i=1}^{k} w_i = 1$和校准分数si，加权分位数可以正式定义为：

$$\hat{q} = \inf\left\{q : \sum_{i=1}^{n} w_i \mathbb{I}\{s_i \lt q\} \geq 1 - \alpha\right\}$$

其中inf是下确界，$\mathbb{I}$是指示函数，条件为真时值为1，否则为0。在此上下文中，下确界是使得不等式成立的最小q值。这只是在前面提到的笔记本中看到的加权分位数的更严格定义。其余过程与之前完全相同。

实际上，我们可以通过几种不同的方式将其用于时间序列问题。例如：

  1. 我们可以考虑一个长度为 K 的滑动窗口，并有一个固定长度的权重向量，长度为 K。在这种方案下，我们将对时间序列中的每个点应用权重，作用于最后 K 个点，并计算该点的预测区间。这些权重可以是等权重，甚至可以是衰减权重，从而捕捉时间因素。
  2. 当我们将多个时间序列一起建模时，我们可以确保权重反映另一个时间序列与我们正在生成区间的时间序列的接近程度，同时考虑时间上下文。

关键在于，我们在设计权重时可以尽可能发挥创意。指导原则是，权重应反映校准数据点与您正在生成预测区间的数据点之间的差异程度。请记住，我们之前讨论过上界，$d_{TV}(s(z), s(z^i))$。我们选择的权重将抵消这一项。当我们对距离我们关心的数据点“较远”的数据点赋予较小权重时，这会降低覆盖差距的上界，使其更加紧凑。

现在，让我们了解另一种非常简单的修改，以应对分布偏移。

自适应共形推断 (ACI)

在2021年，Gibbs等人（参考文献 13）提出了另一种在在线场景下处理分布偏移（特别是时间序列）的方法。时间序列数据通常一次只出现一个数据点，而这种处理分布偏移的方法依赖于这种在线特性，因为它建议根据不断流入的数据持续调整预测区间，从而使预测区间适应变化的分布。这种方法称为 自适应共形推断（ACI），可以与任何预测算法集成，在非平稳条件下提供稳健的预测集。

在传统的共形预测中，我们有一个评分函数，$s(x, y)$，和一个分位数函数，$\hat{q}$，它给出了预测区间，$\hat{C}$。注意，底层不确定性模型（被共形化为 $\hat{q}$）可以是任何估计不确定性的方法，如分位数回归、PDF、MC Dropouts 等。当数据可交换时，我们在校准数据上计算的 $\hat{q}$ 将在未来的测试数据点上保持有效。但当分布发生偏移时，这个 $\hat{q}$ 也会开始变得不太相关。为了解决这个问题，作者建议定期重新估计这些函数，使其与最新的数据观测保持一致。具体来说，在每个时间点 t，基于最新数据拟合一个新的评分函数 $s_t(.)$ 和一个新的分位数函数 $\tilde{q}_t$(.)。

为此，他们定义了 覆盖不足率，$M_t(\alpha)$，作为真实标签 $Y_t$ 落在预测区间 $\hat{C}$ 之外的概率，其中概率是在校准数据和测试数据点上计算的。我们希望 覆盖不足率，$M_t(\alpha)$，等于 $\alpha$（期望误差率）。但由于数据分布正在变化，$M_t(\alpha)$ 预计不会随时间保持不变，并且可能不等于目标水平 $\alpha$。作者假设，对于每个时间点 $t$，可能存在一个最优覆盖水平 $\alpha_t^*$，使得覆盖不足率 $M_t(\alpha_t^*)$ 近似等于 $\alpha$。

为了估计这个 $\alpha_t^*$，作者提出一个简单的在线更新方程。该更新考虑了先前观测的经验覆盖不足率，然后减少或增加我们对 $\alpha_t^*$ 的估计。具体来说，如果我们设置 $\alpha_1 = \alpha$，我们可以将误差定义为：

$$err_t = \begin{cases} 1, & \text{if } Y_t \notin \hat{C}(a_t), \\ 0, & \text{otherwise}, \end{cases}$$

现在，我们可以递归地定义更新步骤为：

$$\alpha_{t+1} = \alpha_t + \gamma(\alpha - err_t)$$

这里，$err_t$ 作为历史误覆盖率估计，$\gamma \gt 0$ 是步长（超参数，后续详述）。因此，当 $err_t = 0$（预测落在区间内）时，$\alpha - err_t$ 为正，进而将 $\alpha_{t+1}$ 更新为高于 $\alpha_t$ 的值。这反过来使得预测区间变窄（根据我们定义的 $\gamma$）。同理，当 $err_t = 1$（预测落在区间外）时，预测区间变宽。

一种自然的替代更新方式，它更多考虑历史信息，是使用过去时间步的加权平均：

$$\alpha_{t+1} = \alpha_t + \gamma \left( \alpha - \sum_{i=1}^{t} w_i \cdot err_i \right)$$

其中 $\{w_i\}_{1 \leq i \leq t} \in [0,1]$ 是递增权重序列，满足 $\sum_{i=1}^{t} w_i = 1$。该更新不是只看最后一个时间步的误覆盖率估计，而是考察近期历史。这使其在理论上至少更稳健一些。论文报告两种策略没有显著差异。他们用于决定权重的一种策略是：

$$w_i = \frac{0.95^{t-s}}{\sum_{s=1}^{t} 0.95^{t-s}}$$

他们报告，通过简单更新和加权更新获得的预测区间轨迹几乎相同，但加权更新的轨迹明显更平滑，$\alpha_t$ 的局部变化更小。

现在，我们花些时间理解步长参数 $\gamma$ 的影响。其直觉与深度学习中的学习率非常相似。$\gamma$ 决定了我们更新 $\alpha$ 的幅度。值越大，更新越快，反之亦然。论文还给出直觉：分布漂移越大，$\gamma$ 的值越大。在他们所有的实验中，他们使用 $\gamma = 0.005$，理由是他们发现这个值能使轨迹相对平滑，同时仍足够大以使 $\alpha_t$ 适应分布漂移。我们可以将其视为控制“适应”强度的参数，使我们能够在非自适应区间和强自适应区间之间的频谱中移动。

现在，让我们看看如何在实际中应用这些。

使用共形预测进行预测

我们没有找到任何现成的实现，涵盖我们想在此展示的所有技术，尤其是具备以下特性的实现：

因此，我们提供了一个文件（src/conformal/conformal_predictions.py），其中包含了与 neuralforecast 预测兼容的必要实现，并拥有统一的 API。它足够简单易懂。我们将讲解代码的主要部分，但若要了解所有部分如何协同工作，你只需查看该文件即可。

我们讨论的所有方法，如回归共形预测、共形分位数回归和共形化不确定性估计，都已使用相同的 API 编码。让我们看看最基本的回归共形预测以理解此 API。它可以在文件中的 ConformalPrediction 类中找到。其余技术都继承此类并稍作修改。所有类的编码方式使得它们接收来自 neuralforecast（或 statsforecast）的预测数据框，并使用相同的命名约定来共形化这些预测。理论上，任何可以转换为预期格式的预测都可以使用这些类。预期格式中的列如下：

除了这些列，我们还会有一列（或多列）预测值，按相应名称命名。

在开始生成共形预测之前，我们还需要一些数据和预测。我们使用本章一直使用的相同数据（M4），并额外创建了一个拆分（校准）。使用新的训练数据，我们使用 level = 90 创建了以下三个预测：

  1. LSTM 点预测（LSTM）
  2. 带分位数回归的 LSTM（LSTM_QR）
  3. 基于PDF（正态分布）的LSTM (LSTM_PDF))

笔记本提醒

要跟随完整代码，请使用Chapter17文件夹中的名为09-Conformal_Techniques.ipynb的笔记本。

笔记本包含完整的代码，但我们可以从已经将数据分割成Y_train_df、Y_calib_df和Y_test_df，并在字典中生成和存储预测的prediction_dict。让我们看看准备好的数据框的前五行，以了解我们正在处理的数据类型。

图17.16：我们正在使用的Y_calib_df的前五行。这是我们编写的共形预测类所期望的格式

现在，让我们言归正传，开始创建预测区间。

回归的共形预测

ConformalPrediction类提供了一种结构化方法，基于所选模型在标定数据集上的预测来计算预测区间。它包括以下输入参数：

使用该类的主要函数有：

这些类是外部 API。内部有一些实际定义实现方式的方法。让我们在常规共形预测的背景下查看主要方法。

calculate_scores 方法定义如下，我们仅使用校准数据集计算绝对残差作为分数：

def calculate_scores(self, Y_calib_df):
    Y_calib_df = Y_calib_df.copy()
    Y_calib_df["calib_scores"] = np.abs(Y_calib_df["y"] - Y_calib_df[self.model])
    return Y_calib_df

get_quantile 方法使用定义的 $\alpha$ 为每个 unique_id 计算分位数：

def get_quantile(self, Y_calib_df):
    def get_qhat(Y_calib_df):
        n_cal = len(Y_calib_df)
        q_level = np.ceil((n_cal + 1) * (1 - self.alpha)) / n_cal
        return np.quantile(
                Y_calib_df["calib_scores"].values, q_level, method="higher"
            )
    return Y_calib_df.groupby(
"unique_id").apply(get_qhat).to_dict()

calc_prediction_interval 方法使用计算出的 q_hat 和均值预测来生成预测区间。对于常规共形预测，过程如下：

def calc_prediction_interval(self, Y_test_df, q_hat):
    return (
            Y_test_df[self.model] - Y_test_df["unique_id"].map(q_hat),
            Y_test_df[self.model] + Y_test_df["unique_id"].map(q_hat),
        )

现在，我们将其用于预测。

from src.conformal.conformal_predictions import ConformalPrediction
Y_calib_df, Y_test_df = prediction_dict['LSTM']
# Y_calib_df & Y_test_df have forecasts in column named "LSTM"
cp = ConformalPrediction(model="LSTM", level=level)
# Calibrating the model
cp.fit(Y_calib_df=Y_calib_df)
# Generating Prediction intervals
Y_test_df_cp = cp.predict(Y_test_df=Y_test_df)

生成的带有预测区间的 DataFrame（图17.15）将有两列——分别为下预测区间和上预测区间的 LSTM-CP-lo-90 和 LSTM-CP-hi-90。CP 是分配给共形预测类的方法标签。

我们可以按如下方式检查任何对象的方法名称：

>> cp.method_name
'Vanilla Conformal Prediction (CP)

我们来看看预测数据框的样子（CP 是分配给共形预测类的方法标签。我们可以通过做 cp.method_name 来检查任何对象的方法名称）：

图17.17：带有预测区间的生成数据框

现在，我们还需要计算区间的覆盖率和平均长度，以评估这些预测区间的效果。我们使用之前相同的方法来实现。先不讨论每种方法的性能，我们将讨论留到最后，现在先着眼于创建预测区间。

那么，我们继续介绍下一项技术。

共形分位数回归

应用这种共形预测方法的第一个条件是，底层分位数回归应已生成一组预测区间。因此，我们使用LSTM_QR中训练的模型。

普通共形预测与CQR的主要区别在于分数计算方式和预测区间。因此，我们可以继承ConformalPrediction，只需重新定义这两个方法。

我们来看一下calculate_scores方法。

def calculate_scores(self, Y_calib_df):
    Y_calib_df = Y_calib_df.copy()
    lower_bounds = Y_calib_df[self.lower_quantile_model]
    upper_bounds = Y_calib_df[self.upper_quantile_model]
    Y_calib_df["calib_scores"] = np.maximum(
            lower_bounds - Y_calib_df["y"], Y_calib_df["y"] - upper_bounds
        )
    return Y_calib_df

我们刚刚实现了前面看到的公式。self.lower_quantile_model和self.upper_quantile_model是CQR已生成区间的列名。

现在，我们还需要定义calc_prediction_interval方法。

def calc_prediction_interval(self, Y_test_df, q_hat):
    return (
            Y_test_df[self.lower_quantile_model] - Y_test_df["unique_id"].map(q_hat),
            Y_test_df[self.upper_quantile_model] + Y_test_df["unique_id"].map(q_hat),
        )

q_hat是一个字典，包含为每个unique_id计算的分位数。因此，我们只需获取CQR现有的预测区间，并通过输入数据帧中的unique_id映射分位数来调整它。

现在，我们用这个进行预测。API与之前完全相同。

from src.conformal.conformal_predictions import ConformalizedQuantileRegression
Y_calib_df, Y_test_df = prediction_dict['LSTM_QR']
# Forecast in column "LSTM_QR"
cp = ConformalizedQuantileRegression(model="LSTM_QR", level=level)
cp.fit(Y_calib_df=Y_calib_df)
Y_test_df_cqr = cp.predict(Y_test_df=Y_test_df)

现在，我们来看讨论过的第三项技术。

共形不确定性估计

如果还记得之前的讨论，要使用这项技术，我们需要一个可以进一步共形化的不确定性估计。这就是我们选择前面用过的另一项技术PDF的原因。但用MC Dropout也可以轻松实现。我们只需要标准偏差或类似能捕捉每个数据点不确定性的指标。

我们为此使用之前生成的 LSTM_PDF 预测。尽管模型预测正态分布的均值和标准差，但其内部用于生成预测区间。因此，我们之前定义的 PDF 模型的输出将是预测区间，但我们想要标准差。别担心。我们知道预测区间是使用正态分布创建的。因此，从预测区间重新推导标准差并不困难。

$$Lower\ Bound = \mu - Z \cdot \sigma$$

$$Upper\ Bound = \mu + Z \cdot \sigma$$

使用基本数学，我们可以推导出：

$$\sigma = \frac{Upper \ Bound - Lower \ Bound}{Z}$$

而 Z 非常容易获得。我们可以使用 scipy.stats.norm 来实现。下面是一种从预测区间获取标准差的方法（请记住，这仅适用于使用正态分布创建的 PDF）。

from scipy.stats import norm
def calculate_standard_deviation(upper_bound, point_prediction, confidence_level):
    # Calculate the Z-value from the confidence level
    z_value = norm.ppf((1 + confidence_level) / 2)
    # Calculate the standard deviation
    sigma = (upper_bound - point_prediction) / z_value
    return sigma
def reverse_engineer_sd(X, model_tag, level):
    X["std"] = calculate_standard_deviation(
        X[f"{model_tag}-hi-{level}"], X[model_tag], level / 100
    )
    return X

现在，我们将其添加到我们的 Y_calib_df 和 Y_test_df 中。

Y_calib_df = reverse_engineer_sd(Y_calib_df, "LSTM_Normal", level)
Y_test_df = reverse_engineer_sd(Y_test_df, "LSTM_Normal", level)

现在，让我们看看如何定义该类。我们需要在此处添加之前不需要的额外信息——不确定性估计的列名。因此，我们定义新类（仍然继承 ConformalPrediction）如下：

class ConformalizedUncertaintyEstimates(ConformalPrediction):
    def __init__(
        self,
        model: str,
        uncertainty_model: str,
        level: Optional[float] = None,
        alias: str = None,
    ):
        super().__init__(model, level, alias)
        self.method = "Conformalized Uncertainty Intervals"
        self._mthd = "CUE"
        self.uncertainty_model = uncertainty_model

我们定义了一个额外的参数 uncertainty_model，并将其他参数传递给父类。

现在，这非常简单。我们需要定义如何计算分数：

def calculate_scores(self, Y_calib_df):
    Y_calib_df = Y_calib_df.copy()
        uncertainty = Y_calib_df[self.uncertainty_model]
    Y_calib_df["calib_scores"] = (
            np.abs(Y_calib_df["y"] - Y_calib_df[self.model]) / uncertainty
        )
    return Y_calib_df

以及 calc_prediction_interval 方法：

def calc_prediction_interval(self, Y_test_df, q_hat:
    uncertainty = Y_test_df[self.uncertainty_model]
    return (
            Y_test_df[self.model] - uncertainty * Y_test_df["unique_id"].map(q_hat),
            Y_test_df[self.model] + uncertainty * Y_test_df["unique_id"].map(q_hat),
        )

就是这样。现在，我们有了一个新类，可以对不确定性估计进行共形化。让我们用它来获取我们一直在处理的数据集的预测。

from src.conformal.conformal_predictions import ConformalizedUncertaintyEstimates
# We have saved uncertainty estimates in "std"
cp = ConformalizedUncertaintyEstimates(model="LSTM_Normal", uncertainty_model="std", level=level)
cp.fit(Y_calib_df=Y_calib_df)
Y_test_df_pdf = cp.predict(Y_test_df=Y_test_df)

现在，我们再来看看之前介绍的两种技术，它们更适合存在分布偏移的时间序列问题。

加权共形预测

我们之前看到，加权共形预测就是在校准数据上应用合适的权重，使得与测试点相似的样本获得比不相似样本更大的权重。关键区别仅在于分位数的计算方式。

这意味着我们可以在底层使用任何共形预测技术，但需要计算加权分位数而非简单分位数。因此，从实现角度来看，我们可以将这个类视为其他已定义技术的包装器，并将它们转换为加权共形预测。

尽管加权共形预测可以通过多种方式实现，接受不同类型的权重（如时间上的权重、unique_id上的权重等），但我们将实现一种更简单的基于回看窗口的加权共形预测。我们选择最近的 K 个时间步，并使用给定的权重在这些 K 个分数上计算加权分位数。权重可以是简单的均匀权重，对所有 K 个步长赋予相同权重，也可以是衰减权重，对最近的分数赋予最高权重，甚至可以是完全自定义的权重。

因此，让我们如下定义类的 __init__：

class WeightedConformalPredictor:
    def __init__(
        self,
        conformal_predictor: ConformalPrediction,
        K: int,
        weight_strategy: str,
        custom_weights: list = None,
        decay_factor: float = 0.5,
    ):
    …

这里，K 是窗口，conformal_predictor 是底层使用的共形预测类（应为我们已定义的三个类之一）。我们可以将权重策略定义为 uniform、decay 或 custom，分别对应均匀权重、衰减权重或自定义权重。decay_factor 决定衰减权重策略中权重衰减的速度，custom_weights 允许你精确指定这些 K 个时间步上的权重。

虽然我们不会在正文中展示全部代码，但会过一遍关键部分，以便你理解发生了什么。但我强烈建议你花些时间消化文件中的代码。

首先，是我们的 fit 方法。在这个方法中，我们仅使用底层共形预测器的分数计算并存储校准数据框。

def fit(self, Y_calib_df):
    self.calib_df = self.conformal_predictor.calculate_scores(
            Y_calib_df.sort_values(["unique_id", "ds"])
        )

现在，让我们看看 predict 方法的主要部分。

def predict(self, Y_test_df):
    # Groupby unique_id
    …
    # Calculate quantiles for each unique_id
    self.q_hat = {}
    for unique_id, group in grouped_calib:
        # Take the last K timesteps
        group = group.iloc[-self.K :]
        scores = group["calib_scores"].values
        # Calculate weights based on the last K timesteps
        total_timesteps = len(scores)
        weights = self._calculate_weight(total_timesteps)
        normalized_weights = weights / weights.sum()
        # Calculate quantile for the current unique_id
        quantile = self.get_weighted_quantile(
                scores, normalized_weights, self.conformal_predictor.alpha
            )
        self.q_hat[unique_id] = quantile
    # Calculate prediction intervals using the underlying conformal predictor's method
    lo, hi = self.conformal_predictor.get_prediction_interval_names()
    Y_test_df[lo], Y_test_df[hi] = (
 self.conformal_predictor.calc_prediction_interval(Y_test_df, self.q_hat)
        )
    return Y_test_df

现在，取我们之前看到的一个方法，在其上应用加权共形预测包装器。在我们的示例中，选择简单的 ConformalPrediction。让我们看看如何使用这个类：

from src.conformal.conformal_predictions import WeightedConformalPredictor
Y_calib_df, Y_test_df = prediction_dict['LSTM']
# Defining an underlying conformal predictor
cp = ConformalPrediction(model="LSTM", level=level)
# using the defined conformal predictor in weighted version
weighted_cp = WeightedConformalPredictor(
    conformal_predictor=cp,
    K=50,
    weight_strategy="uniform",
)
weighted_cp.fit(Y_calib_df=Y_calib_df)
Y_test_df_wcp = weighted_cp.predict(Y_test_df=Y_test_df)

这将创建带有标签 CP_Wtd 的预测区间。我们可以通过执行 weighted_cp.method_name 来检查标签。

现在，该实现有一个小缺点。尽管我们在分数中考虑了时间顺序，但我们仍然有一个固定的校准集。因此，在利用最新数据点“重新拟合”或校准之前，我们仍然在使用相同的校准数据集。所以，仔细想来，这个方法理想情况下应该以在线方式应用，即每次预测新的时间步时，前一个时间步（带有实际值）应该被添加到校准数据中。我们还提供了一个能够以在线方式完成此操作的替代实现。我们不深入实现细节，因为核心逻辑相同，但API不同，使得更新校准数据成为可能。完整实现可在文件中的OnlineWeightedConformalPredictor类中找到。

让我们看看如何使用它。首先，我们定义设置，初始化类，并拟合校准。

from src.conformal.conformal_predictions import OnlineWeightedConformalPredictor
cp = ConformalPrediction(model="LSTM", level=level)
online_weighted_cp = OnlineWeightedConformalPredictor(
    conformal_predictor=cp,
    K=50,
    weight_strategy="uniform",
)
online_weighted_cp.fit(Y_calib_df=Y_calib_df)
joblib.dump(online_weighted_cp, "path/to/saved/file.pkl")

现在，在推理过程中，我们可以对每个时间步执行类似操作：

# Loading the saved model
online_weighted_cp = joblib.load("path/to/saved/file.pkl")
# current timestep data = current
# past timestep actuals = last_timestep_actuals
prediction = online_weighted_cp.predict_one(current_test)
# updating the calibration data using the last timestep actuals
online_weighted_cp.update(last_timestep_actuals)

对于我们在知道实际值的测试数据上进行评估的特殊情况，还有另一种方法可以对数据进行类似的在线预测：offline_predict。

Y_test_df_wcpo = online_weighted_cp.offline_predict(Y_test_df=Y_test_df)

现在，让我们看最后一个方法。

自适应共形推断

最后，我们来看自适应共形推断。这也可以作为其他共形预测方法的包装器实现，因为该技术涉及更新$\alpha$，以便在分布偏移的情况下保持覆盖范围。由于该技术的性质，我们只能以在线方式应用它，即利用可用数据在每个时间步更新alpha。因此，它将具有与我们之前看到的OnlineWeightedConformalPredictor相同的API。

完整的类可在src/conformal/conformal_predictions.py中找到，但这里我们将查看一些主要部分以便理解。让我们先看看__init__函数：

class OnlineAdaptiveConformalInference:
    def __init__(
        self,
        conformal_predictor: ConformalPrediction,
        gamma: float = 0.005,
        update_method: str = "simple",
        momentum_bw: float = 0.95,
        per_unique_id: bool = True,
    ):
        …

类似于WeightedConformalPredictor，我们接受一个底层共形预测器(conformal_predictor)。此外，我们还有gamma，即步长($\gamma$)；update_method，可以是simple（仅使用最后一个时间步进行更新）或momentum（使用误差轨迹的移动平均）。最后，我们还可以定义momentum_bw，即用于计算过去误差加权平均的动量回权重。较高的动量（例如0.95）使轨迹更平滑，使过去的误覆盖率衰减更慢。另一个极端（0.05）使得加权平均更灵敏，更接近“简单”方法。最后，我们还有一个参数，用于对每个unique_id分别计算误差或将所有误差合并。

像往常一样，我们有一个fit方法，它使用校准数据集计算分数并使其就绪。我们也可以将其视为在线实现中的预热期。$\alpha$更新将使用校准数据作为初始历史记录来启动。

def fit(self, Y_calib_df):
    """
    Fit the conformal predictor model with calibration data.
    """
    self.calib_df = self.conformal_predictor.calculate_scores(Y_calib_df)
    self.scores_by_id = (
        self.calib_df.groupby("unique_id")["calib_scores"].apply(list).to_dict()
    )
    # Some more code to initialize necessary data structures
    …
    return self

现在，让我们看看predict_one方法，它预测一个时间步。

def predict_one(self, current_test):
    unique_ids = current_test["unique_id"].unique()
    predictions = []
    for unique_id in unique_ids:
        group_scores = self.scores_by_id.get(unique_id, [])
        if group_scores:
            # Determine the appropriate alpha to use
            alpha = (
                self.alphat[unique_id]
                if self.per_unique_id
                else self.alphat_global
            )
            # Calculate quantile for the current unique_id
            self.q_hat = {
                unique_id: np.quantile(group_scores, 1 - alpha, method="higher")
            }
            # Calculating prediction intervals using conformal_predictor
            …
            current_test[lo] = lower.values
            current_test[hi] = upper.values
            # Storing most recent prediction
            …
        # Collecting and returning concatenated predictions
        …

代码注释充分，你应该能理解。对于每个unique_id，我们获取分数，获取适当的$\alpha$，计算分位数，利用这些信息使用底层共形预测器计算预测区间，并存储预测结果以备后用。

现在，我们来看一下 update 方法，该方法可以在时间步的实际值可用时使用。

def update(self, new_data):
    new_scores = self.conformal_predictor.calculate_scores(new_data)
    for unique_id, score in zip(
        new_scores["unique_id"], new_scores["calib_scores"]
    ):
        # Updating score trajectory with new score
    …
        # Retrieve stored predictions and calculate adapt_err
        if self.per_unique_id:
            lower, upper = self.predictions[unique_id]
            actual_y = new_data.loc[new_data["unique_id"] == unique_id, "y"].values[0]
            adapt_err = int(actual_y < lower or actual_y > upper)
            # Update alpha updates the alpha using simple or momentum method
            self.update_alpha(unique_id, adapt_err)
        else:
            # Do the same update at a global error-pooled way

让我们看看如何使用它。首先，我们定义设置、初始化类并拟合校准。

from src.conformal.conformal_predictions import OnlineAdaptiveConformalInference
cp = ConformalPrediction(model="LSTM", level=level)
aci_cp = OnlineAdaptiveConformalInference(
    conformal_predictor=cp,
    gamma=0.005,
    update_method="simple",
)
aci_cp.fit(Y_calib_df=Y_calib_df)
joblib.dump(aci_cp, "path/to/saved/file.pkl")

现在，在推理过程中，我们可以对每个时间步执行类似以下操作：

# Loading the saved model
aci_cp = joblib.load("path/to/saved/file.pkl")
# current timestep data = current
# past timestep actuals = last_timestep_actuals
prediction = aci_cp.predict_one(current_test)
# updating the calibration data using the last timestep actuals
aci_cp.update(last_timestep_actuals)

对于我们在测试数据上进行评估的特定情况（我们知道实际值），还有另一种方法可以对数据进行类似的在线预测：offline_predict。

Y_test_df_aci = aci_cp.offline_predict(Y_test_df=Y_test_df)

既然我们已经看到了所有技术的实际应用，让我们来看看它们的表现如何。

评估结果

如果你一直跟随笔记本进行操作，你会知道我们为每种技术计算了覆盖率（coverage）和平均长度（average length）。回顾一下，覆盖率衡量的是实际值落在预测区间内的时间百分比，而平均长度衡量的是预测区间的平均宽度。对于 level=90，我们期望覆盖率大约为 90% 或 0.9，并且平均长度尽可能小，同时不损害覆盖率。

让我们看看以下总结覆盖率和平均长度的图表：

图 17.18：所有共形技术的覆盖率（按 unique_id 分），颜色表示它们接近 0.9 的程度（我们设置了 level=90）

图 17.19：所有共形技术的平均长度（按 unique_id 分），颜色表示它们的大小

以下是图例的快速回顾，以理解不同的列，这些列只是根据应用组合了以下标签：

开门见山，我们可以看到，针对分布偏移进行校正的技术表现最佳。两个表格右侧（接近0.9覆盖率和较低平均长度）有更多“绿色”。记住，我们从常规共形预测开始，然后针对分布偏移进行了校正。我们可以注意到，对于大多数 unique_id，常规共形预测对 level = 90 的区间过宽。覆盖率大于0.9，多数情况下为1.0。但加权校正（CP_Wtd）和自适应共形推断（ACI）都使预测区间更窄，且更接近我们的预期水平。两者之间没有明确的赢家，这需要根据你的数据集来评估。

至此，我们关于概率预测的讨论就结束了。希望借此你能自信地进入这个鲜为人知的领域，并为你的预测对象带来更多价值。

现在，让我们从非常高的层面审视时间序列预测中一些很少受到关注、但在许多领域却相当相关的细分主题。

时间序列预测中的冷门路径

秉承罗伯特·弗罗斯特《未选择的路》的精神，本节探讨时间序列预测中鲜为人知但极具影响力的技术。正如弗罗斯特选择了人迹罕至的那条路，我们深入研究那些虽非主流、却在各领域提供独特见解和潜在突破的专业方法。

间歇性或稀疏需求预测

间歇性时间序列预测处理的是零星、不规则事件的数据，常见于零售业，其中某些产品可能销售不频繁。它对于管理库存、避免缺货或过度库存至关重要，尤其是慢周转商品。传统方法难以应对这些模式，因为对于大多数此类物品，需求期望趋于零，因此需要专门技术来提高准确性，这使其成为零售预测的必备工具。

在此，我们将快速列举一些专为间歇性预测设计的替代算法及其实现位置。

可解释性

我们在第10章《全局预测模型》中已经介绍了一些机器学习模型的可解释性技术。虽然其中一些技术（如SHAP和LIME）仍可应用于深度学习模型，但没有一种在设计中考虑时间因素。这是因为所有这些技术都是为更通用的目的（如分类和回归）开发的。尽管如此，在深度学习模型和时间序列模型的可解释性方面已有一些工作。在此，我将列出几篇有前景的、直面时间因素的论文：

虽然这不是一个详尽的列表，但这些是我们认为重要且有前景的一些工作。这是一个活跃的研究领域，随着时间推移，新技术将会出现。

冷启动预测

冷启动预测解决了预测无历史数据产品需求的挑战，这是零售、制造和消费品等行业中的常见问题。它在推出新产品、引入品牌或扩展到新地区时出现。传统的统计预测模型如 ARIMA 或指数平滑无法处理这一问题，因为它们需要大量历史数据。

但并非毫无希望。我们确实有一些方法可以处理此类情况：

分层预测

分层预测处理的是可以分解为嵌套层级的时间序列，例如产品类别或地理区域。这些结构要求预测保持一致性，即较低层级的预测应汇总到较高层级，反映聚合关系。分组时间序列结合了不同的层次结构（例如产品类型和地理位置），增加了复杂性。目标是产生与数据自然聚合一致且准确的预测，这使得分层预测对于管理跨多个维度的大量时间序列的企业至关重要。

有一些技术可以对所有预测进行分解、聚合或协调，使它们以合理的方式相加。Rob Hyndman 的免费时间序列预测圣经（预测：原理与实践）中的第10章对此进行了详细讨论，是快速入门的好资源（温馨提示：数学内容较多）。你可以在这里找到该章节：https://otexts.com/fpp2/hierarchical.html。对于更实际的方法，你可以查看 Nixtla 的 hierarchicalforecast 库：https://nixtlaverse.nixtla.io/hierarchicalforecast/index.html。

至此，本书最长的章节之一告一段落。恭喜你坚持下来并消化了所有信息。请随意使用笔记本，尝试不同的代码和选项，以更好地理解其中的内容。

小结

在本章中，我们探讨了多种生成概率预测的技术，如概率密度函数、分位数函数、蒙特卡洛丢弃法和共形预测。我们深入研究了每种方法，并学习了如何将其应用于真实数据。对于共形预测这一活跃研究领域，我们学习了不同的方式来通过底层机制（如共形化分位数回归、共形化不确定性估计等）对预测区间进行共形校正。最后，我们介绍了一些使共形预测在时间序列问题上更有效的技巧。

最后，我们还涉及了一些冷门主题，如间歇性需求预测、可解释性、冷启动预测和层次预测。

在本书的下一部分，我们将探讨一些预测机制，如多步预测、交叉验证和评估。

参考文献

以下是本章的参考文献：

  1. Tony Duan, A. Avati, D. Ding, S. Basu, A. Ng, and Alejandro Schuler. (2019). NGBoost: Natural Gradient Boosting for Probabilistic Prediction. International Conference on Machine Learning. https://proceedings.mlr.press/v119/duan20a/duan20a.pdf..
  2. Y. Gal and Zoubin Ghahramani. (2015). Dropout as a Bayesian Approximation: Representing Model Uncertainty in Deep Learning. International Conference on Machine Learning. https://proceedings.mlr.press/v48/gal16.html..
  3. Valentin Flunkert, David Salinas, and Jan Gasthaus. (2017). DeepAR: Probabilistic Forecasting with Autoregressive Recurrent Networks. International Journal of Forecasting. https://www.sciencedirect.com/science/article/pii/S0169207019301888..
  4. Koenker, Roger. (2005). Quantile Regression. Cambridge University Press. pp. 146–7. ISBN 978-0-521-60827-5. http://www.econ.uiuc.edu/~roger/research/rq/QRJEP.pdf..
  5. Spyros Makridakis, Evangelos Spiliotis, Vassilios Assimakopoulos. (2020). The M4 Competition: 100,000 time series and 61 forecasting methods. International Journal of Forecasting. https://www.sciencedirect.com/science/article/pii/S0169207019301128..
  6. Loic Le Folgoc 和 Vasileios Baltatzis 和 Sujal Desai 和 Anand Devaraj 和 Sam Ellis 和 Octavio E. Martinez Manzanera 和 Arjun Nair 和 Huaqi Qiu 和 Julia Schnabel 和 Ben Glocker。（2021）。Is MC Dropout Bayesian?。arXiv preprint arXiv: Arxiv-2110.04286。https://arxiv.org/abs/2110.04286。
  7. Nicolas Dewolf 和 Bernard De Baets 和 Willem Waegeman。（2023）。Valid prediction intervals for regression problems。Artif. Intell. Rev.https://doi.org/10.1007/s10462-022-10178-5。
  8. V. Vovk, A. Gammerman, 和 G. Shafer。（2005）。Algorithmic learning in a random world。Springerhttps://link.springer.com/book/10.1007/b106715。
  9. Anastasios N. Angelopoulos 和 Stephen Bates。（2021）。A Gentle Introduction to Conformal Prediction and Distribution-Free Uncertainty Quantification。arXiv preprint arXiv: Arxiv-2107.07511。https://arxiv.org/abs/2107.07511。
  10. Harris Papadopoulos, Kostas Proedrou, Volodya Vovk, 和 Alex Gammerman。（2002）。Inductive confidence machines for regression。Machine Learning: ECML 2002。https://doi.org/10.1007/3-540-36755-1_29。
  11. Romano, Yaniv 和 Patterson, Evan 和 Candes, Emmanuel。（2019）。Conformalized Quantile Regression。Advances in Neural Information Processing Systems。https://proceedings.neurips.cc/paper_files/paper/2019/file/5103c3584b063c431bd1268e9b5e76fb-Paper.pdf。
  12. R. Barber, E. Candès, Aaditya Ramdas, 和 R. Tibshirani。（2022）。Conformal prediction beyond exchangeability。Annals of Statistics。https://projecteuclid.org/journals/annals-of-statistics/volume-51/issue-2/Conformal-prediction-beyond-exchangeability/10.1214/23-AOS2276.full。
  13. Gibbs, Isaac 和 Candes, Emmanuel。（2021）。Adaptive Conformal Inference Under Distribution Shift。Advances in Neural Information Processing Systems。https://proceedings.neurips.cc/paper_files/paper/2021/file/0d441de75945e5acbc865406fc9a2559-Paper.pdf。

延伸阅读

要了解更多关于本章涵盖的主题，请查看以下资源。

留下评论！

感谢您从Packt Publishing购买本书——希望您喜欢！您的反馈非常宝贵，有助于我们改进和成长。阅读完成后，请花一点时间在亚马逊上留下评论；只需一分钟，但对像您这样的读者意义重大。

扫描二维码或访问链接以免费获取您选择的一本电子书。

https://packt.link/NzOWQ

第4部分

预测机制

在最后这一部分中，我们将介绍几个对于构建行业级预测系统至关重要的概念。我们讨论一些鲜有提及的概念，如多步预测，并深入探讨评估预测的细节。

本部分包含以下章节：