Python 机器学习算法交易实用指南（三）

原文：zh.annas-archive.org/md5/fcb09c483bdb21866eb6782158d1f8d5

译者：飞龙

协议：CC BY-NC-SA 4.0

第八章：时间序列模型

在上一章中，我们重点关注了适用于横截面数据的线性模型，其中输入数据属于与其旨在解释或预测的输出相同的时间段。在本章中，我们将关注时间序列数据，其中观察值按周期而不同，这也创建了自然的排序。我们的目标是识别数据中的历史模式，并利用这些模式来预测时间序列在未来的行为。

在上一章中，我们已经遇到了既有横截面又有时间序列维度的面板数据，并了解了法马-麦克白斯回归如何估计随时间和跨资产采取某些因子风险的价值。然而，随时间的回报之间的关系通常相当低，因此这个过程可能会大部分忽视时间维度。本章的模型侧重于时间序列模型，其中过去的价值包含关于未来发展的预测信号。时间序列模型还可以预测随后用于横截面模型的特征。

具体来说，在本章中，我们将重点关注从先前观察到的数据中提取信号以预测同一时间序列未来值的模型。交易的时间维度使得将时间序列模型应用于市场、基本和替代数据变得非常流行。随着越来越广泛的连接设备收集可能包含预测信号的定期测量，时间序列数据将变得更加普遍。关键应用包括资产回报、相关性或协方差、或波动率的预测。

在本章中，我们将重点关注线性时间序列模型，作为本书第四部分应用于时间序列数据的非线性模型（如递归或卷积神经网络）的基线。我们首先介绍诊断时间序列特征的工具，包括平稳性，并提取捕捉潜在模式的特征。然后，我们介绍单变量和多变量时间序列模型，并将其应用于宏观数据和波动性模式的预测。最后，我们介绍协整概念以及如何应用它来开发配对交易策略。

具体来说，我们将涵盖以下主题：

如何使用时间序列分析来诊断通知建模过程的诊断统计数据
如何估计和诊断自回归和移动平均时间序列模型
如何构建自回归条件异方差（ARCH）模型来预测波动性
如何构建向量自回归模型
如何利用协整性进行配对交易策略

诊断和特征提取的分析工具

时间序列数据是一系列由离散时间间隔分隔的值，通常是均匀间隔的（除了缺失值）。

两个不同时间点 t[i]、t[j]之间的周期数Δt= t[i] - t[j]称为滞后，对于每个时间序列有 T-1 个滞后。

如果一个时间序列是独立同分布的随机变量ε[t]的序列，并且具有有限的均值和方差，则称其为白噪声。

如果一个时间序列可以写成过去扰动ε[t]的加权和，这些扰动也称为创新，这里假设它们代表白噪声，以及序列的均值μ，则该时间序列是线性的。

时间序列分析的一个关键目标是理解由系数 a[i]驱动的动态行为。

在本章的大多数示例中，我们使用由美联储提供的数据，您可以使用我们在第二章中介绍的pandas datareader进行访问，市场和基本数据。本节的代码示例可在笔记本tsa_and_arima中找到。

如何分解时间序列模式

时间序列数据通常包含各种模式的混合，可以分解为几个组成部分，每个部分代表一种基本模式类别。特别是，时间序列通常由系统组件趋势、季节性和周期性以及非系统噪声组成。这些组件可以在加法线性模型中组合，特别是当波动不依赖于系列水平时，或者在非线性乘法模型中组合。

这些组件可以自动拆分。statsmodels包括一种简单的方法，将时间序列拆分为趋势、季节和残差分量，使用移动平均。我们可以将其应用于工业制造生产的月度数据，该数据具有强烈的趋势和季节性分量，如下所示：

import statsmodels.tsa.api as tsa
industrial_production = web.DataReader('IPGMFN', 'fred', '1988', '2017-12').squeeze()
components = tsa.seasonal_decompose(industrial_production, model='additive')
ts = (industrial_production.to_frame('Original')
      .assign(Trend=components.trend)
      .assign(Seasonality=components.seasonal)
      .assign(Residual=components.resid))
ts.plot(subplots=True, figsize=(14, 8));

结果图显示了加法组件。假设趋势和季节性组件更确定且更适合简单外推，则残差组件将成为额外建模的焦点：

还有更复杂的基于模型的方法，这些方法包含在 GitHub 上可用的参考资料中。

如何计算滚动窗口统计量

考虑到时间序列数据的顺序排列，自然而然地计算给定长度周期的熟悉描述性统计量，以检测稳定性或行为变化，并获得捕获系统方面的平滑表示，同时过滤噪声。

滚动窗口统计为此过程服务：它们生成一个新的时间序列，其中每个数据点表示原始数据的一定期间内计算的摘要统计量。移动平均是最熟悉的例子。原始数据点可以以相等的权重进入计算，或者使用权重来，例如，强调最近的数据点。指数移动平均递归地计算权重，使得远离现在的数据点权重缩小或衰减。新的数据点通常是所有前面数据点的摘要，但也可以从周围的窗口计算。

pandas库包括非常灵活的功能来定义各种窗口类型，包括滚动窗口、指数加权窗口和扩展窗口。在第二步中，您可以对由窗口捕获的每个数据应用计算。这些计算包括针对单个系列的内置标准计算，例如平均值或总和，对于几个系列的相关性或协方差，以及用户定义的函数。以下部分中的移动平均和指数平滑示例使用了这些工具。

移动平均和指数平滑

早期的预测模型包括带有指数权重的移动平均模型，称为指数平滑模型。我们将在后面的线性时间序列中再次遇到移动平均作为关键的构建块。

依赖指数平滑方法的预测使用过去观测的加权平均值，其中随着观测时间的推移，权重呈指数衰减。因此，更近期的观测获得更高的关联权重。这些方法在时间序列不具有非常复杂或突然的模式时很受欢迎。

指数平滑是一种流行的技术，基于过去观测的加权平均，随着观测时间的推移，权重呈指数衰减。换句话说，观测越近，相关的权重越高。这个框架可以快速生成可靠的预测，适用于广泛的时间序列，这是一个巨大的优势，对于工业应用至关重要。

如何测量自相关

自相关（也称为串行相关）将相关概念调整为时间序列环境：正如相关系数衡量两个变量之间线性关系的强度一样，自相关系数ρ[k]衡量相隔给定滞后 k 的时间序列值之间的线性关系的程度：

因此，我们可以计算时间序列中 T-1 个滞后的每个自相关系数；T 是序列的长度。自相关函数（ACF）根据滞后计算相关系数。

对于大于 1 的滞后自相关（即，超过一个时间步的观测之间的），反映了这些观测之间的直接相关以及中间数据点的间接影响。偏自相关消除了这种影响，只测量给定滞后距离的数据点之间的线性依赖关系。**偏自相关函数（PACF）**提供了一旦较短滞后的相关性效应被移除后产生的所有相关性。

有算法根据 PACF 与 ACF 之间的确切理论关系来估计样本自相关的偏自相关。

自相关图简单地是顺序滞后 k=0,1,...,n 的 ACF 或 PACF 的图。它允许我们一眼就检查各个滞后的相关结构。自相关图的主要用途是在去除确定性趋势或季节性影响后检测任何自相关。ACF 和 PACF 都是线性时间序列模型设计的关键诊断工具，我们将在时间序列转换的以下部分中审查 ACF 和 PACF 图的示例。

如何诊断并实现平稳性

稳定时间序列的统计特性，如均值、方差或自相关，与时期无关，即它们随时间不变化。因此，稳态意味着时间序列没有趋势或季节效应，而且当为不同的滚动窗口计算描述统计量（如均值或标准差）时，它们是常数或随时间变化不大的。它回归到其均值，偏差具有恒定的幅度，而短期波动在统计意义上总是相同的。

更正式地说，严格稳态要求任何时间序列观测子集的联合分布与时间无关，而与所有时刻的矩一致。因此，除了均值和方差外，高阶矩，如偏度和峰度，也需要保持不变，而不考虑不同观测之间的滞后。在大多数应用中，我们将稳态限制为一阶和二阶矩，以使时间序列具有恒定的均值、方差和自相关。

请注意，我们特别允许不同滞后观测之间的依赖关系，就像我们希望线性回归的输入数据与结果相关联一样。稳态意味着这些关系是稳定的，这有助于预测，因为模型可以专注于学习发生在稳定统计属性内的系统模式。这很重要，因为经典统计模型假设时间序列输入数据是稳态的。

以下部分介绍了帮助检测数据非稳态的诊断方法，以及帮助满足这些假设的转换方法。

时间序列转换

要满足线性时间序列模型的稳态假设，我们需要转换原始时间序列，通常需要几个步骤。常见的转换包括将（自然）对数应用于将指数增长模式转换为线性趋势并稳定方差。通货紧缩意味着将一个时间序列除以另一个导致趋势行为的系列，例如将名义系列除以价格指数以将其转换为实际测量。

如果一个系列恢复到一个稳定的长期线性趋势，则它是趋势稳定的。通常可以通过使用线性回归拟合趋势线并使用残差，或者通过将时间指数包括为回归或 AR(I)MA 模型的独立变量（请参阅下一节有关单变量时间序列模型的内容）来使其稳定，可能结合对数化或通货紧缩。

在许多情况下，去趋势化不足以使序列平稳。相反，我们需要将原始数据转换为周期间和/或季节间差异的系列。换句话说，我们使用相邻数据点或季节滞后值相互减去的结果。请注意，当将这种差分应用于对数变换的序列时，结果代表金融背景下的瞬时增长率或收益率。

如果单变量序列在差分 d 次后变为平稳，则称其为 d 阶整合，如果 d=1，则称其为简单整合。这种行为是由所谓的单位根引起的。

如何诊断和解决单位根问题

单位根对确定将时间序列变为平稳的转换提出了特殊问题。时间序列通常被建模为我们将更详细地探讨的具有以下自回归形式的随机过程，作为 ARIMA 模型的构建块：

其中当前值是过去值的加权和加上随机扰动。这样的过程具有以下形式的特征方程：

如果此方程的根之一等于 1，则该过程被称为具有单位根。它将是非平稳的，但不一定需要趋势。如果特征方程的剩余根的绝对值小于 1，则过程的一阶差分将是平稳的，并且过程是整合的（阶数为 1）或 I(1)。如果附加根的绝对值大于 1，则积分的阶数较高，并且将需要额外的差分。

在实践中，利率或资产价格的时间序列通常不是平稳的，例如，因为不存在价格水平使序列回归。非平稳序列最突出的例子是给定起始价格 p[0]（例如，股票的首次公开发行价格）和满足以下条件的白噪声扰动 ε 的价格时间序列的随机游走：

重复代换表明，当前值 p[t] 是所有先前扰动或创新 ε 以及初始价格 p[0] 的总和。如果方程包括一个常数项，则随机游走被称为具有漂移。因此，随机游走是以下形式的自回归随机过程：

具有特征方程，该方程具有单位根并且既不稳定也不是一阶积分。一方面，鉴于 ε 的 i.i.d. 性质，时间序列的方差等于 σ²，这不是二阶稳定的，并且暗示着原则上，该系列随时间可能采用任何变量。另一方面，进行第一次差分，Δp[t]=p[t]-p[t-1]，留下 Δp[t]=ε [t]，这是稳定的，鉴于对 ε 的统计假设。

单位根非平稳系列的定义特征是长期记忆：因为当前值是过去扰动的总和，所以大的创新持续的时间比均值回归、平稳的系列长得多。

除了使用相邻数据点之间的差异来消除恒定变化模式之外，还可以应用季节性差分来消除季节性变化的模式。这涉及到以表示季节模式长度的滞后距离为单位，相隔四个季度或 12 个月的值之间的差异，以消除季节性和线性趋势。

识别正确的转换，尤其是差分的适当数量和滞后期并不总是明确的。一些规则已经被提出，总结如下：

10+滞后期的正自相关: 可能需要更高阶的差分。
滞后期接近零或为负，或一般较小且无规律: 不需要更高阶的差分。
滞后期-1 的自相关< -0.5: 系列可能存在过多的差分。
轻微的过度或不足的差分可以通过 AR 或 MA 项进行修正。
最优差分通常会产生最低的标准差，但并非总是如此。
无差分的模型假设原始系列是平稳的，包括均值回归。通常包括一个常数项，以允许非零均值。
具有一次差分次序的模型假设原始系列具有恒定的平均趋势，并且应包括一个常数项。
具有两个差分次序的模型假设原始系列具有时变趋势，不应包括一个常数。

一些作者推荐分数差分作为使集成系列保持稳定的更灵活的方法，可能能够在离散间隔中保留比简单或季节性差异更多的信息或信号（参见 GitHub 上的参考资料）。

单位根检验

统计单位根检验是一种常见的客观确定（额外）差分是否必要的方法。这些是用于确定是否需要差分的平稳性的统计假设检验。

增广 Dickey-Fuller (ADF) 检验评估一个时间序列样本是否具有单位根的零假设，对比平稳性的备择假设。它对时间趋势、第一滞后和所有滞后差异的差异化时间序列进行回归，并根据滞后时间序列值的系数值计算一个检验统计量。 statsmodels 使得实现变得容易（参见伴随的笔记本）。

形式上，时间序列 y[t] 的 ADF 测试运行线性回归：

其中 α 是一个常数，β 是时间趋势的系数，p 是模型中使用的滞后数。α=β =0 约束意味着一个随机游走，而只有 β=0 意味着有漂移的随机游走。滞后顺序通常使用在第七章中介绍的 AIC 和 BIC 信息标准来决定，线性模型。

ADF 检验统计量使用样本系数 γ，在单位根非平稳的零假设下等于零，在其他情况下为负。它旨在证明，对于一个积分序列，滞后的系列值不应该在预测超出滞后差异的第一差异时提供有用的信息。

如何应用时间序列变换

以下图表显示了纳斯达克股票指数和工业生产在原始形式下的 30 年时间序列，以及应用对数和随后应用第一和季节性差异（在滞后 12 ）之后的转换版本。图表还显示了 ADF p 值，这使得我们可以在两种情况下拒绝单位根非平稳的假设：

我们可以进一步分析转换后时间序列的相关特征，使用一个 Q-Q 图来比较时间序列观测值的分布的分位数和正态分布的分位数，以及基于 ACF 和 PACF 的自相关图。

对于纳斯达克图表，我们注意到虽然没有趋势，但方差不是恒定的，而是在 1980 年代末、2001 年和 2008 年的市场动荡时期附近显示出集群的尖峰。Q-Q 图突出显示了分布的尾部肥大，极端值比正态分布所暗示的更频繁。ACF 和 PACF 显示出类似的模式，多个滞后处的自相关显著：

对于工业生产的月度时间序列，我们注意到 2008 年危机后出现了一个大的负异常值，以及 Q-Q 图中相应的偏斜。自相关远高于纳斯达克回报，并且平稳下降。PACF 在滞后 1 和 13 处显示出明显的正自相关模式，并且在滞后 3 和 4 处有显著的负系数：

单变量时间序列模型

多元线性回归模型将感兴趣的变量表达为预测变量或输入变量的线性组合。单变量时间序列模型将感兴趣时间点的时间序列值与系列滞后值的线性组合以及可能的过去干扰项相关联。

虽然指数平滑模型是基于数据中的趋势和季节性的描述，但 ARIMA 模型旨在描述数据中的自相关性。ARIMA(p, d, q) 模型需要稳定性，并利用两个构建块：

自回归 (AR) 项由时间序列的 p 滞后值组成
移动平均 (MA) 项包含 q 滞后干扰

I 代表积分，因为模型可以通过对系列进行 d 次微分来考虑单位根非平稳性。自回归术语强调 ARIMA 模型意味着时间序列对其自身值的回归。

我们将介绍 ARIMA 构建块，简单的自回归 (AR) 和移动平均 (MA) 模型，并解释如何将它们结合在自回归移动平均 (ARMA) 模型中，该模型可以作为 ARIMA 模型考虑系列集成，或者包含外生变量作为 AR(I)MAX 模型。此外，我们将说明如何包含季节性 AR 和 MA 项以扩展工具箱，以便还包括 SARMAX 模型。

如何构建自回归模型

顺序 p 的 AR 模型旨在捕捉不同滞后时间序列值之间的线性相关性，并可写为：

这与 y[t] 的滞后值的多重线性回归非常相似。该模型具有以下特征方程：

解此方程的 x 的解的倒数是特征根，如果这些根在绝对值上都小于 1，则 AR(p) 过程是稳定的，否则是不稳定的。对于稳态系列，多步预测将收敛到系列的平均值。

我们可以使用熟悉的最小二乘法估计模型参数，使用 p+1，..., T 观测数据以确保每个滞后项和结果都有数据。

如何确定滞后数量

实际上，挑选适当的滞后项顺序 p 是一个挑战。用于串行相关性的时间序列分析工具起着关键作用。ACF 估计不同滞后观测之间的自相关性，这反过来是直接和间接线性相关性的结果。

因此，对于阶数为 k 的 AR 模型，ACF 将显示出一直到滞后 k 的显著串行相关，并且由于线性关系间接影响所导致的惯性，将延伸到后续滞后，最终因效应减弱而逐渐消失。另一方面，PACF 仅度量观察值之间的直接线性关系，间隔给定滞后，因此不会反映超过 k 的滞后的相关性。

如何诊断模型拟合度

如果模型捕获了滞后之间的线性依赖关系，则残差应该类似于白噪声。

除了检查 ACF 以验证自相关系数的显著性缺失外，Ljung-Box Q 统计量还允许我们测试残差序列是否遵循白噪声的假设。零假设是所有 m 个串行相关系数都为零，而备择假设是一些系数不为零。测试统计量是从不同滞后 k 的样本自相关系数 ρ[k] 计算得到，并且遵循一个 Χ² 分布：

正如我们将看到的，statsmodels 提供了关于不同滞后的系数显著性的信息，不显著的系数应予以删除。如果 Q 统计量拒绝了无自相关的原假设，则应考虑额外的 AR 项。

如何构建移动平均模型

阶数为 q 的 MA 模型使用过去 q 个扰动，而不是时间序列的滞后值，如下：

由于我们没有观察到白噪声扰动值 ε[t]，MA(q) 不像我们迄今为止见过的那样是一个回归模型。MA(q) 模型不是使用最小二乘法估计的，而是使用最大似然估计（MLE），或者在系列开始时初始化或估计扰动，然后递归和迭代地计算剩余部分。

MA(q)模型的名称源于将每个 y[t] 的值表示为过去 q 个创新的加权移动平均。换句话说，当前的估计值代表了相对于模型过去错误的修正。与指数平滑或估计季节性时间序列分量不同，MA(q)模型中移动平均的使用旨在预测未来值，而不是去噪声或估计过去值的趋势周期。

MA(q) 过程始终是平稳的，因为它们是白噪声变量的加权和，而这些变量本身是平稳的。

如何确定滞后数

由 MA(q) 过程产生的时间序列由前 q 个模型预测的残差驱动。因此，MA(q) 过程的 ACF 将显示出到滞后 q 的值的显著系数，然后急剧下降，因为这是假定系列值生成的方式。

AR 和 MA 模型之间的关系

AR(p) 模型可以使用重复替代表示为 MA(∞) 过程。当对其系数的大小施加约束时，一个 MA(q) 过程，它变得可逆，并且可以表示为 AR(∞) 过程。

如何构建 ARIMA 模型及其扩展

自回归积分移动平均 ARIMA(p, d, q) 模型结合了 AR(p) 和 MA(q) 过程，利用这些构建模块的互补性，并通过使用更紧凑的形式简化模型开发，并减少参数的数量，从而减少过拟合的风险。

该模型还通过使用时间序列值的 d^(th) 差分来消除单位根非平稳性。ARIMA(p, 1, q) 模型与使用系列的第一次差分的 ARMA(p, q) 模型相同。使用 y' 表示非季节性差分 d 次后的原始系列，ARIMA(p, d, q) 模型就是：

ARIMA 模型也是使用最大似然估计的。根据实现的不同，高阶模型通常会包含低阶模型。例如，statsmodels 包括所有低阶 p 和 q 项，并且不允许移除小于最高值的滞后项系数。在这种情况下，高阶模型总是拟合得更好。小心不要通过使用过多的项过度拟合模型。

如何确定 AR 和 MA 项的数量

由于 AR(p) 和 MA(q) 项相互作用，ACF 和 PACF 提供的信息不再可靠，只能作为起点使用。

传统上，AIC 和 BIC 信息准则一直被用来在选择模型设计时依赖于样本内拟合。或者，我们可以依赖于样本外测试来交叉验证多个参数选择。

以下摘要提供了在考虑 AR 和 MA 模型时选择模型阶数的一些通用指导：

超过 PACF 截断的滞后是所示的 AR 项的数量。如果差分序列的 PACF 显示出明显的截断和/或滞后 1 的自相关是正的，则添加一个或多个 AR 项。
超过 ACF 截断的滞后是所示的 MA 项的数量。如果差分序列的 ACF 显示出明显的截断和/或滞后 1 的自相关是负的，则考虑向模型中添加一个 MA 项。
AR 和 MA 项可能会相互抵消彼此的影响，因此如果您的模型同时包含两者，尝试总是将 AR 和 MA 项的数量减少 1，以避免过度拟合，特别是如果更复杂的模型需要超过 10 次迭代才能收敛。
如果 AR 系数之和接近于 1，并且表明模型的 AR 部分存在单位根，则应消除 1 个 AR 项并再次进行差分。
如果 MA 系数之和接近于 1，并且表明模型的 MA 部分存在单位根，则应消除 1 个 MA 项，并将差分阶数降低 1。
不稳定的长期预测表明模型的 AR 或 MA 部分可能存在单位根。

添加特征 - ARMAX

ARMAX 模型在 ARMA 时间序列模型的右侧添加输入变量或协变量（假设序列是平稳的，因此我们可以跳过差分）：

这类似于线性回归模型，但很难解释，因为 β 对 y[t] 的影响不是像线性回归中 x[t] 增加一个单位那样简单。相反，在方程的右侧存在 y[t] 的滞后值意味着系数只能在考虑到响应变量的滞后值时进行解释，这几乎是不直观的。

添加季节性差分 - SARIMAX

对于具有季节效应的时间序列，我们可以包括捕捉季节性周期性的 AR 和 MA 项。例如，当使用月度数据且季节效应长度为一年时，季节性 AR 和 MA 项将反映这个特定的滞后长度。

然后 ARIMAX(p, d, q) 模型变成了 SARIMAX(p, d, q) x (P, D, Q)[s] 模型，这样写起来稍微复杂一些，但是 GitHub 上的参考资料，包括 statsmodels 的文档，提供了详细信息。

现在我们将使用宏观数据构建一个季节性 ARMA 模型来说明实施过程。

如何预测宏观基本面

我们将为 1988-2017 年间的工业生产时间序列构建一个 SARIMAX 模型。正如在关于分析工具的第一部分中所说明的，数据已经进行了对数变换，我们正在使用季节性（滞后-12）差分。我们使用一个包括常规 AR 和 MA 参数范围的滚动窗口对 10 年的训练数据进行模型估计，并评估 1 步预测的 RMSE，如下所示的简化代码（详情请参阅 GitHub）：

for p1 in range(4):                # AR order
    for q1 in range(4):            # MA order
        for p2 in range(3):        # seasonal AR order
            for q2 in range(3):    # seasonal MA order
                y_pred = []
                for i, T in enumerate(range(train_size, len(data))):
                    train_set = data.iloc[T - train_size:T]
                    model = tsa.SARIMAX(endog=train_set,            # model specification
                                        order=(p1, 0, q1),
                                        seasonal_order=(p2, 0, q2, 12)).fit()

                    preds.iloc[i, 1] = model.forecast(steps=1)[0]    # 1-step ahead forecast

                mse = mean_squared_error(preds.y_true, preds.y_pred)
                test_results[(p1, q1, p2, q2)] = [np.sqrt(mse),
                                                  preds.y_true.sub(preds.y_pred).std(),
                                                  np.mean(aic)]

我们还收集了 AIC 和 BIC 标准，显示了 0.94 的非常高的秩相关系数，其中 BIC 倾向于比 AIC 稍少参数的模型。RMSE 最佳的五个模型为：

                 RMSE         AIC         BIC
p1 q1 p2 q2                                  
2  3  1  0   0.009323 -772.247023 -752.734581
3  2  1  0   0.009467 -768.844028 -749.331586
2  2  1  0   0.009540 -770.904835 -754.179884
   3  0  0   0.009773 -760.248885 -743.523935
   2  0  0   0.009986 -758.775827 -744.838368

我们重新估计了一个 SARIMA(2, 0, 3) x (1, 0, 0) 模型，如下所示：

best_model = tsa.SARIMAX(endog=industrial_production_log_diff, order=(2, 0, 3),
                         seasonal_order=(1, 0, 0, 12)).fit()
print(best_model.summary())

我们得到以下摘要：

系数是显著的，并且 Q 统计量拒绝了进一步自相关的假设。相关图表同样表明我们已成功消除了序列的自相关性：

如何使用时间序列模型来预测波动性

单变量时间序列模型的一个特别重要的应用领域是波动率的预测。金融时间序列的波动率通常随时间而变化，但是波动率的变化会聚集在一起。方差的变化对于使用经典的 ARIMA 模型进行时间序列预测构成挑战。为了解决这一挑战，我们现在将对波动率进行建模，以便我们可以预测方差的变化。

异方差性是变量方差变化的技术术语。自回归条件异方差（ARCH）模型将误差项的方差表达为前期误差的函数。更具体地说，它假定误差方差遵循一个 AR(p) 模型。

广义自回归条件异方差（GARCH）模型将范围扩展到 ARMA 模型。时间序列预测通常将 ARIMA 模型用于时间序列的预期均值，并将 ARCH/GARCH 模型用于时间序列的预期方差。2003 年诺贝尔经济学奖授予了罗伯特·恩格尔和克莱夫·格兰杰，以表彰他们开发了这类模型。前者还在纽约大学斯特恩商学院（参见 GitHub 引用）运营着波动性实验室，提供了许多关于我们将讨论的模型及其众多扩展的在线示例和工具。

自回归条件异方差（ARCH）模型

ARCH(p) 模型简单地是将 AR(p) 模型应用于时间序列模型的残差方差，该模型使得时间 t 的方差条件于方差的滞后观测。更具体地说，误差项 ε[t] 是原始时间序列上线性模型（如 ARIMA）的残差，并分为一个时间相关的标准差 σ[t] 和一个扰动 z[t]，如下所示：

ARCH(p) 模型可以使用 OLS 进行估计。Engle 提出了一种使用拉格朗日乘子检验来确定适当的 ARCH 阶数的方法，该检验对应于线性回归中所有系数为零的假设的 F 检验（参见第七章，线性模型）。

该模型的一个优点是它产生波动率，并估计正的超额峰度——即，相对于正态分布有厚尾，这与关于回报的经验观察一致。弱点包括模型假设正和负波动冲击具有相同的效应，因为它依赖于前一冲击的平方，而资产价格已知对正和负冲击有不同的反应。ARCH 模型也没有提供关于金融时间序列变化源的新见解，因为它只是机械地描述条件方差。最后，ARCH 模型很可能会过度预测波动性，因为它对回报系列的大规模孤立冲击反应缓慢。

对于正确规范的 ARCH 模型，标准化残差（除以标准偏差期间的模型估计）应类似于白噪声，并且可以进行 Ljung-Box Q 检验。

泛化 ARCH —— GARCH 模型

ARCH 模型相对简单，但通常需要许多参数来捕获资产收益系列的波动性模式。广义 ARCH（GARCH）模型适用于对数收益系列 r[t]，具有扰动项 ε[t] = r[t ]- μ，如果遵循 GARCH(p, q) 模型：

GARCH(p, q) 模型假设误差项 ε[t] 的方差服从 ARMA(p, q) 模型。

类似于 ARCH 模型，GARCH(1,1) 过程的尾部分布比正态分布更重。该模型面临与 ARCH 模型相同的弱点。例如，它对正负冲击的反应一样。

选择滞后阶数

要为 ARCH 和 GARCH 模型配置滞后阶数，使用训练以预测原始系列均值的时间序列的平方残差。残差是零中心化的，因此它们的平方也是方差。然后检查平方残差的 ACF 和 PACF 图以识别时间序列方差中的自相关模式。

如何构建波动率预测模型

为资产收益序列开发波动率模型包括四个步骤：

基于 ACF 和 PACF 显示的序列依赖性构建金融时间序列的 ARMA 时间序列模型。
再次依靠残差序列的 ACF 和 PACF 检查模型的残差是否具有 ARCH/GARCH 效应。
如果序列相关效应显著，则指定一个波动率模型，并联合估计均值和波动率方程。
仔细检查拟合的模型，如有必要，进行精细调整。

当将波动率预测应用于回报序列时，序列依赖性可能有限，因此可以使用恒定均值代替 ARMA 模型。

arch 库提供了几个选项来估计波动率预测模型。它提供了几个选项来建模预期均值，包括一个常数均值，上面关于单变量时间序列模型的 AR(p) 模型以及更近期的异质自回归过程（HAR），它使用每日（1 天）、每周（5 天）和每月（22 天）的滞后来捕捉短期、中期和长期投资者的交易频率。

可以与几个条件异方差模型一起定义和估计均值模型，除了 ARCH 和 GARCH 外，还包括允许正负收益之间的非对称效应的指数 GARCH（EGARCH）模型，以及异质 ARCH（HARCH）模型，它补充了 HAR 均值模型。

我们将使用 1998 年至 2017 年的每日纳斯达克收益来演示 GARCH 模型的使用情况（有关详细信息，请参阅笔记本 arch_garch_models）：

nasdaq = web.DataReader('NASDAQCOM', 'fred', '1998', '2017-12-31').squeeze()
nasdaq_returns = np.log(nasdaq).diff().dropna().mul(100) # rescale to facilitate optimization

重新缩放的日收益率系列仅具有有限的自相关性，但与均值的平方偏差具有显着的记忆，反映在缓慢衰减的 ACF 和 PACF 上，前两者的值很高，并且仅在前六个滞后之后才截断：

plot_correlogram(nasdaq_returns.sub(nasdaq_returns.mean()).pow(2), lags=120, title='NASDAQ Daily Volatility')

函数plot_correlogram生成了以下输出：

因此，我们可以估计一个 GARCH 模型来捕捉过去波动率的线性关系。我们将使用滚动 10 年窗口来估计一个 GARCH(p, q)模型，其中 p 和 q 的范围为 1-4，以生成 1 步的外样本预测。然后，我们将比较预测波动性的 RMSE 与实际收益偏离其均值的平方之间的 RMSE，以确定最具预测性的模型。我们使用修剪数据来限制极端收益值的影响，这些值反映在波动性的非常高的正偏度中：

trainsize = 10 * 252  # 10 years
data = nasdaq_returns.clip(lower=nasdaq_returns.quantile(.05),
                           upper=nasdaq_returns.quantile(.95))
T = len(nasdaq_returns)
test_results = {}
for p in range(1, 5):
    for q in range(1, 5):
        print(f'{p} | {q}')
        result = []
        for s, t in enumerate(range(trainsize, T-1)):
            train_set = data.iloc[s: t]
            test_set = data.iloc[t+1]  # 1-step ahead forecast
            model = arch_model(y=train_set, p=p, q=q).fit(disp='off')
            forecast = model.forecast(horizon=1)
            mu = forecast.mean.iloc[-1, 0]
            var = forecast.variance.iloc[-1, 0]
            result.append([(test_set-mu)**2, var])
        df = pd.DataFrame(result, columns=['y_true', 'y_pred'])
        test_results[(p, q)] = np.sqrt(mean_squared_error(df.y_true, df.y_pred))

GARCH(2, 2)模型实现了最低的 RMSE（与 GARCH(4, 2)的值相同，但参数较少），因此我们继续估计此模型以检查摘要：

am = ConstantMean(nasdaq_returns.clip(lower=nasdaq_returns.quantile(.05),
                                      upper=nasdaq_returns.quantile(.95)))
am.volatility = GARCH(2, 0, 2)
am.distribution = Normal()
model = am.fit(update_freq=5)
print(model.summary())

输出显示了最大化的对数似然以及通常在选择基于样本内表现的模型时要最小化的 AIC 和 BIC 准则（见第七章，线性模型）。它还显示了均值模型的结果，本例中只是一个常数估计，以及常数 omega 的 GARCH 参数、AR 参数α和 MA 参数β，所有这些参数都具有统计学意义：

现在让我们探讨多个时间序列模型和协整概念，这将实现新的交易策略。

多元时间序列模型

多元时间序列模型旨在同时捕捉多个时间序列的动态，并利用这些序列之间的依赖关系进行更可靠的预测。

方程组

像我们刚刚讨论的 ARMA 方法这样的单变量时间序列模型仅限于目标变量与其滞后值或滞后扰动以及在 ARMAX 情况下的外生序列之间的统计关系。相比之下，多元时间序列模型还允许其他时间序列的滞后值影响目标。这种效应适用于所有系列，导致复杂的相互作用，如下图所示：

除了可能更好的预测之外，多变量时间序列还用于获得跨系列依赖关系的洞察。例如，在经济学中，多变量时间序列用于了解对某一变量的政策变化，例如利率，可能在不同的视角下如何影响其他变量。我们将要看的多变量模型产生的冲击-响应函数就是为此目的而服务的，并且允许我们模拟一个变量如何对其他变量的突然变化做出响应。格兰杰因果分析的概念分析了一个变量是否对另一个变量的预测有用（在最小二乘意义上）。此外，多变量时间序列模型允许对预测误差方差进行分解，以分析其他系列的贡献。

向量自回归 (VAR) 模型

我们将看到向量自回归 VAR(p) 模型如何通过创建包含所有 k 系列的 p 拖后值的 k 方程系统来扩展 AR(p) 模型。在最简单的情况下，k=2 的 VAR(1) 模型如下所示：

这个模型可以用矩阵形式更加简洁地表示：

自身滞后的系数提供了有关系列动态的信息，而交叉变量系数则提供了有关系列间交互的一些见解。此符号扩展到 k 系列和阶数 p，如下所示：

VAR(p) 模型还需要平稳性，以便从单变量时间序列建模的初始步骤延续下来。首先，探索系列并确定必要的转换，然后应用增广迪基-富勒检验来验证每个系列是否满足平稳性标准，并在否则应用进一步的转换。它可以通过有条件的 OLS 进行估计，初始信息或最大似然，对于正态分布的误差是等效的，但对于其他情况不是。

如果一些或所有的 k 系列是单位根非平稳的，则它们可能是协整的。将单位根概念扩展到多个时间序列意味着两个或多个系列的线性组合是平稳的，因此是均值回归的。VAR 模型没有能力处理这种情况而不进行差分，而应使用向量误差校正模型 (VECM，请参阅 GitHub 上的参考资料)。我们将进一步探讨协整性，因为如果存在并且被假定为持续存在，它可以用于配对交易策略。

滞后阶数的确定也从每个系列的 ACF 和 PACF 中获得线索，但受到同一滞后阶数适用于所有系列的约束。在模型估计之后，残差诊断也要求产生类似白噪声的结果，并且模型选择可以使用样本内信息准则或更好地使用样本外预测性能来交叉验证备选模型设计，如果最终目标是使用模型进行预测。

如在单变量情况中提到的那样，对原始时间序列的预测要求我们在训练模型之前撤消应用于使系列稳态的转换。

如何使用 VAR 模型进行宏观基本面预测

我们将扩展单变量示例，将一系列工业生产的月度数据和美联储数据服务提供的一系列消费者情绪的月度时间序列加入，我们将使用熟悉的pandas-datareader库从 1970 年到 2017 年检索数据：

df = web.DataReader(['UMCSENT', 'IPGMFN'], 'fred', '1970', '2017-12').dropna()
df.columns = ['sentiment', 'ip']

对工业生产系列进行对数转换，并使用这两个系列的滞后 12 个月进行季节性差分产生稳态结果：

df_transformed = pd.DataFrame({'ip': np.log(df.ip).diff(12),
                              'sentiment': df.sentiment.diff(12)}).dropna()

test_unit_root(df_transformed) # see notebook for details and additional plots

          p-value
ip          0.0003
sentiment   0.0000

这使我们得到以下系列：

为了限制输出大小，我们将仅使用前 480 个观测值使用statsmodels的VARMAX实现（允许使用可选的外生变量）来估计一个 VAR(1)模型，使用恒定的趋势：

model = VARMAX(df_transformed.iloc[:480], order=(1,1), trend='c').fit(maxiter=1000)

这导致以下摘要：

输出包含两个时间序列方程的系数，如前述的 VAR(1)示例所述。statsmodels 提供了诊断图来检查残差是否满足白噪声假设，在这种简单情况下并不完全符合：

外样本预测可按以下方式生成：

preds = model.predict(start=480, end=len(df_transformed)-1)

实际值和预测值的可视化显示了预测滞后于实际值，并且无法很好地捕捉非线性的外样本模式：

共整合 - 具有共同趋势的时间序列

综合多变量系列的概念因该过程的所有组成系列可能是单独集成的，但该过程在一个或多个线性组合产生新的稳态系列的意义上并非联合集成而复杂化。

换句话说，两个共整合系列的组合具有一个稳定的均值，这个线性组合会恢复到该均值。具有这种特性的多变量系列被称为共整合。当个别系列积分的次序较高且线性组合减少了积分的总次序时，这也适用。

协整与相关性不同：两个系列可能高度相关，但不一定是协整的。例如，如果两个增长系列是彼此的恒定倍数，它们的相关性将很高，但任何线性组合也将增长而不是回归到均值。

仍然可以将 VAR 分析应用于使用 VAR 模型的误差校正形式的集成过程，该模型使用各个系列的一阶差异加上水平的误差校正项。

测试协整

测试协整有两种主要方法：

Engle-Granger 两步法
Johansen 程序

Engle-Granger 方法涉及将一个系列回归到另一个系列，然后对回归残差应用 ADF 单位根检验。如果能拒绝原假设，假设残差是平稳的，那么这些系列就是协整的。这种方法的一个关键优势是，回归系数表示使组合平稳的乘数，即均值回归。当利用协整进行配对交易策略时，我们将回到这个方面。另一方面，这种方法仅限于识别成对系列的协整，而不是更大范围的系列组合。

相比之下，Johansen 程序测试了前一节讨论的矢量自回归（VAR）模型受协整约束的限制。更具体地说，从通用 VAR(p) 前述方程两边减去目标向量后，我们得到误差校正模型（ECM）公式

结果修改后的 VAR(p) 方程只有一个向量项在水平上，即不使用算子 Δ 表示差异。协整的性质取决于系数矩阵 Π 的特性，特别是其秩。虽然这个方程在结构上与 ADF 测试设置相似，但现在存在多个共同趋势和积分阶数的潜在星座，因为涉及到多个系列。有关详细信息，请参阅 GitHub 上列出的参考文献，包括关于个别系列尺度化的实际挑战。

如何利用协整进行配对交易策略

配对交易依赖于两个资产价格之间的稳定、均值回归关系。换句话说，两个价格之间的比率或差异，也称为价差，随时间可能会发散，但最终应返回到相同水平。鉴于这样的一对，该策略包括做多（即购买）表现不佳的资产，因为它需要一段时间的表现优异来填补差距。同时，应做空价格向正方向移动而远离价格锚点的资产，以筹集购买资产的资金。

协整正是由一个共同均值锚定的两个价格系列之间的这种稳定关系。假设协整持续存在，最终必须发生收敛，无论是表现不佳的股票上涨还是表现良好的股票下跌。该策略无论如何都将是有利可图的，这具有对一般市场动态进行对冲的额外优势。

然而，价差将不断变化，有时扩大，有时缩小，或者保持不变，因为两种资产同时移动。配对交易的挑战在于通过调整相对持仓来维持对冲头寸，因为价差变化。

在实践中，鉴于一组资产，配对交易策略将通过对每对运行统计测试来寻找协整对。这里的关键挑战是考虑到多重测试偏差，如第六章，机器学习工作流中概述的那样。statsmodels库实现了 Engle-Granger 协整测试和 Johansen 测试。

为了估算价差，运行线性回归以获得两个集成资产价格系列的线性组合的系数，以产生一个平稳的组合系列。如前所述，使用线性回归来估算系数被称为协整的 Engle-Granger 检验。

总结

在本章中，我们探讨了用于单个系列的单变量情况的线性时间序列模型以及用于几个交互系列的多变量模型。我们遇到了预测宏观基本面的应用，用于风险管理的预测资产或组合波动性的模型以及捕获多个宏观系列动态的多元 VAR 模型，以及协整的概念，这是支撑流行的配对交易策略的基础。

类似于前一章，我们看到线性模型为模型增加了许多结构，即它们做出了强有力的假设，这些假设可能需要转换和广泛测试来验证这些假设是否成立。如果是这样，模型的训练和解释就是直接的，而且模型提供了一个很好的基准案例，更复杂的模型可能会改进，我们将在接下来的章节中看到。

第九章：贝叶斯机器学习

在本章中，我们将介绍贝叶斯方法在机器学习中的应用，以及它们对开发和评估算法交易策略时的不同不确定性视角的增值。

贝叶斯统计学使我们能够量化对未来事件的不确定性，并以原则性的方式在新信息到来时优化我们的估计。这种动态方法很好地适应了金融市场的发展性质。当存在较少相关数据且我们需要系统地整合先验知识或假设时，它特别有用。

我们将看到，贝叶斯方法使得对统计指标、参数估计和预测周围的不确定性有更丰富的见解。应用范围从更精细的风险管理到动态更新的预测模型，其中包含了市场环境的变化。资产配置的黑-利特曼方法（见第五章，《策略评估》）可以解释为贝叶斯模型。它计算预期收益，作为市场均衡和投资者观点的加权平均值，每个资产的波动性，跨资产的相关性以及对每个预测的信心。

具体来说，在本章中，我们将涵盖以下主题：

贝叶斯统计如何应用于机器学习
如何使用 PyMC3 进行概率编程
如何定义和训练机器学习模型
如何运行最先进的抽样方法进行近似推断
如何应用贝叶斯机器学习来计算动态夏普比率，构建贝叶斯分类器和估计随机波动性

本章的参考文献、附加材料链接和代码示例位于 GitHub 存储库相应目录中。请按照第一章提供的安装说明进行操作，《交易的机器学习》。

贝叶斯机器学习的工作原理

经典统计学也被称为频率派，因为它将概率解释为长期内事件的相对频率，即在观察了大量试验之后。在概率的背景下，一个事件是一个实验的一个或多个基本结果的组合，比如两个骰子掷出六个相等的结果中的任何一个，或者某个资产价格在某一天下跌 10%或更多。

贝叶斯统计学相反，将概率视为事件发生的信心或信念的度量。贝叶斯概率的观点为主观观点留下了更多的空间，因此，与频率派解释相比，意见之间的差异更大。这种差异在很少发生的事件中最为显著，以至于无法得出客观的长期频率度量。

换句话说，频率学派统计假设数据是来自人群的随机样本，并旨在识别生成数据的固定参数。相反，贝叶斯统计将数据视为已知的，并认为参数是随机变量，其分布可以从数据中推断出来。因此，频率学派方法要求的数据点至少与要估计的参数一样多。另一方面，贝叶斯方法与较小的数据集兼容，并且非常适合逐个样本进行在线学习。

贝叶斯观点对于许多在某些重要方面罕见或独特的现实事件非常有用。例如，下次选举的结果或市场是否会在三个月内崩溃的问题。在每种情况下，既有相关的历史数据，又有随着事件临近而展开的独特情况。

首先，我们将介绍贝叶斯定理，该定理通过将先验假设与新的经验证据相结合，并将得到的参数估计与频率学派的对应估计进行比较，以晶化通过更新信念来更新概念的过程。然后，我们将演示两种贝叶斯统计推断的方法，这些方法能够揭示潜在参数的后验分布，即未观察到的参数，在不同情况下的预期值等：

共轭先验通过提供闭合形式的解决方案来促进更新过程，但确切的分析方法并不总是可用。
近似推断模拟了由假设和数据组合而成的分布，并使用该分布的样本来计算统计洞察。

如何从经验证据更新假设

牧师托马斯·贝叶斯在 250 多年前提出的定理利用基本的概率理论规定了概率或信念在相关新信息到达时应该如何变化。以下约翰·梅纳德·凯恩斯的引述体现了贝叶斯主义的思维方式：

“当事实发生变化时，我改变我的想法。先生，你会怎么做？”

它依赖于条件概率和全概率以及链式法则；有关这些概念的评论，请参阅 GitHub 上的参考资料。

信念涉及单个或一组参数 θ（也称为假设）。每个参数可以是离散的或连续的。θ 可以是一个一维统计量，比如（离散的）分类变量的模式，或者（连续的）均值，也可以是一个更高维度的值集，比如一个协方差矩阵或深度神经网络的权重。

频率学派统计的一个关键区别在于，贝叶斯假设被表达为概率分布，而不是参数值。因此，虽然频率学派的推断关注点估计，贝叶斯推断则产生概率分布。

贝叶斯定理通过计算从以下输入中得到的后验概率分布来更新对感兴趣参数的信念，如下图所示：

先验分布指示我们考虑每个可能的假设的可能性有多大。
似然函数输出在给定θ参数的某些值的情况下观察到数据集的概率。
证据度量观察到的数据在所有可能的假设下的可能性。因此，它对所有参数值都是相同的，用于将分子标准化：

贝叶斯定理

后验是先验和似然的乘积，除以证据，反映了假设的更新概率分布，同时考虑了先前的假设和数据。从不同的角度看，先验和似然的乘积来自于将数据和参数的联合分布因子分解的链规则的应用。

对于高维、连续变量，制定变得更加复杂，涉及到（多个）积分。一种替代的制定方法使用赔率来表示后验赔率，作为先验赔率乘以似然比的乘积（有关更多细节，请参见参考资料）。

精确推理：最大后验估计

将贝叶斯规则的实际应用于准确计算后验概率的情况非常有限，因为计算分母中的证据项非常具有挑战性。证据反映了在所有可能的参数值上观察到的数据的概率。它也被称为边际似然，因为它需要通过添加或积分参数的分布来对参数的分布进行边际化。这通常只在具有少量假设值的少量离散参数的简单情况下才可能。

**最大后验概率（MAP）**估计利用了证据是一个常数因子，将后验缩放以满足概率分布的要求。由于证据不依赖于θ，后验分布与似然和先验的乘积成比例。因此，MAP 估计选择使后验最大化的θ的值，考虑到观察到的数据和先验信念，即后验的模态。

MAP 方法与最大似然估计（MLE）不同，MLE 定义了概率分布。MLE 选择使观察到的训练数据的似然函数最大化的参数值θ。

从定义的角度看，MAP 与 MLE 的不同之处在于包括了先验分布。换句话说，除非先验是一个常数，否则 MAP 估计θ将与其 MLE 对应物不同：

最大似然估计解往往反映了频率主义的概率估计应该反映观察到的比例的概念。另一方面，先验对 MAP 估计的影响通常相当于将反映先验假设的数据添加到 MLE 中。例如，一个强烈的先验，即硬币有偏的先验可以通过添加偏斜的试验数据来融入 MLE 背景。

先验分布是贝叶斯模型的重要组成部分。我们现在将介绍一些方便的选择，以便进行分析推断。

如何选择先验

先验应反映参数分布的知识，因为它影响 MAP 估计。如果先验不确定，我们需要进行选择，通常从几个合理的选项中选择。一般来说，证明先验的合理性并通过测试替代是否得出相同结论是一个好的做法。

有几种类型的先验：

客观先验最大化数据对后验的影响。如果参数分布未知，我们可以选择一个无信息的先验，比如在参数值的相关范围内称为平坦先验的均匀分布。
相反，主观先验旨在将模型外部的信息纳入估计中。
经验性先验结合了贝叶斯和频率主义方法，利用历史数据消除主观性，例如通过估计各种时刻以适应标准分布。

在机器学习模型的背景下，先验可以被视为一种正则化器，因为它限制了后验可以假设的值。例如，具有零先验概率的参数不是后验分布的一部分。一般来说，更多的好数据可以得出更强的结论并减少先验的影响。

如何保持推断简单 - 共轭先验

当结果后验与先验具有相同类型的分布，只是参数不同时，先验分布与似然的共轭性。当先验和似然都是正态分布时，后验也是正态分布的。

先验和似然的共轭性暗示了后验的闭合形式解，从而便于更新过程并避免使用数值方法来近似后验。此外，由此产生的后验可以用作下一个更新步骤的先验。

让我们使用一个股价波动的二元分类示例来说明这个过程。

如何动态估计资产价格波动的概率

当数据由具有某种成功概率的二元伯努利随机变量组成时，重复试验中的成功次数遵循二项分布。共轭先验是支持区间[0, 1]上的 Beta 分布，并具有两个形状参数，用于对成功概率进行任意先验分布建模。因此，后验分布也是一个 Beta 分布，我们可以通过直接更新参数来得到。

我们将收集不同大小的二元化日度标准普尔 500 指数收益率样本，其中正面结果是价格上涨。从一个不含信息的先验开始，该先验将每个可能的成功概率在区间[0, 1]内分配相等的概率，我们计算不同证据样本的后验概率。

下面的代码示例显示了更新只是简单地将观察到的成功和失败数量添加到先验分布的参数中以获得后验分布：

n_days = [0, 1, 3, 5, 10, 25, 50, 100, 500]
outcomes = sp500_binary.sample(n_days[-1])
p = np.linspace(0, 1, 100)

# uniform (uninformative) prior
a = b = 1
for i, days in enumerate(n_days):
    up = outcomes.iloc[:days].sum()
    down = days - up
    update = stats.beta.pdf(p, a + up , b + down)

以下图表中绘制了结果后验分布。它们说明了从将所有成功概率视为同等可能的均匀先验到越来越尖峰的分布的演变。

经过 500 个样本，概率集中在 2010 年至 2017 年间正面走势的实际概率约为 54.7%。它还显示了 MLE 和 MAP 估计之间的小差异，后者倾向于稍微朝向均匀先验的期望值，如下图所示：

后验概率

在实践中，共轭先验的使用仅限于低维情况。此外，简化的 MAP 方法避免了计算证据项，但即使在其可用时也具有几个缺点；它不返回分布，因此我们无法推导出不确定性的度量，或将其用作先验。因此，我们需要采用数值方法和随机模拟而不是精确推理，我们将在下文介绍。

近似推理：随机与确定性方法

对于大多数实际相关的模型，将无法通过分析方法推导出精确的后验分布并计算潜在参数的期望值。模型可能具有太多的参数，或者后验分布可能对于分析解而言过于复杂。对于连续变量，积分可能没有封闭形式的解，而空间的维数和被积函数的复杂性可能会阻止数值积分。对于离散变量，边缘化涉及对隐藏变量的所有可能配置求和，虽然原则上这总是可能的，但在实践中，我们经常发现可能存在指数多个隐藏状态，因此精确计算是非常昂贵的。

尽管对于某些应用而言，对未观察到的参数的后验分布可能是感兴趣的，但通常主要是要评估期望，例如，进行预测。在这种情况下，我们可以依赖近似推断：

基于马尔可夫链蒙特卡罗（MCMC）抽样的随机技术已经在许多领域推广了贝叶斯方法的使用。它们通常具有收敛到精确结果的能力。在实践中，抽样方法可能需要大量计算，并且通常仅限于小规模问题。
确定性方法，即变分推断或变分贝叶斯，基于对后验分布的解析近似，并且可以很好地扩展到大型应用程序。它们做出简化假设，例如，后验在特定方式上因子化，或者具有特定的参数形式，例如高斯分布。因此，它们不会生成精确结果，并且可以用作抽样方法的补充。

基于抽样的随机推断

抽样是关于从给定分布*p(x)中抽取样本，X=(x[1], ..., x[n])。假设样本是独立的，大数定理确保对于增长的样本数量，给定实例x[i]*在样本中的比例（对于离散情况）对应于其概率，p(x=x[i])。在连续情况下，类似的推理适用于样本空间的给定区域。因此，对样本的平均值可以用作分布参数的期望值的无偏估计。

实际挑战在于确保独立抽样，因为分布是未知的。相关样本可能仍然是无偏的，但倾向于增加估计的方差，因此将需要更多的样本来获得与独立样本一样精确的估计。

从多元分布中抽样在计算上是具有挑战性的，因为随着维数的增加，状态数量呈指数增长。许多算法简化了这个过程（请参阅参考文献以获取概述）。现在，我们将介绍几种基于 MCMC 的流行变体方法。

马尔可夫链蒙特卡罗抽样

马尔可夫链是描述随机漫步的动态随机模型，由转移概率连接的一组状态组成。马尔可夫性质规定该过程没有记忆，并且下一步仅取决于当前状态。换句话说，它在于当前状态的条件下，过去、现在和未来是独立的，即过去状态的信息不会帮助预测未来超出我们从现在所知道的内容。

蒙特卡罗方法依赖于重复的随机抽样来近似可能是确定性的结果，但不允许解析的精确解。它是在曼哈顿计划期间开发的，用于估计原子级别的能量，并获得了其持久的代号以确保保密性。

许多算法将蒙特卡洛方法应用于马尔可夫链，并通常按以下方式进行：

从当前位置开始。
从提议分布中抽取一个新的位置。
在考虑数据和先验分布的情况下评估新位置的概率：
1. 如果足够可能，移动到新的位置。
2. 否则，保持当前位置不变。
从步骤 1 开始重复。
经过一定数量的迭代后，返回所有接受的位置。

MCMC 旨在识别和探索后验的有趣区域，这些区域集中在显著的概率密度上。当它持续移动到后验的附近高概率状态时，无记忆的过程被认为是收敛的，其中接受率增加。一个关键挑战是平衡对样本空间的随机探索的需要和降低接受率的风险。

此过程的初始步骤可能更反映出起始位置而不是后验，并且通常被丢弃为burn-in样本。 MCMC 的一个关键特性是在一定数量的迭代后，过程应该忘记其初始位置。

剩余的样本被称为过程的轨迹。假设收敛，则样本的相对频率近似于后验，可以根据大数定律计算期望值。

正如之前所指出的，估计的精度取决于随机游走收集的样本的串行相关性，每个样本设计上仅取决于前一个状态。更高的相关性限制了对后验的有效探索，并需要经过诊断测试。

设计这样的马尔可夫链的一般技术包括 Gibbs 采样、Metropolis-Hastings 算法和更近期的哈密顿 MCMC 方法，这些方法往往表现更好。

Gibbs 采样

Gibbs 采样将多变量采样简化为一系列一维抽样。从一个起始点开始，它迭代地将 n-1 个变量保持不变，同时抽样第 n 个变量。它将这个样本合并并重复。

该算法非常简单易实现，但产生高度相关的样本，导致收敛速度减慢。其顺序性也阻止了并行化。

Metropolis-Hastings 采样

Metropolis-Hastings 算法基于其当前状态随机提出新的位置，以有效地探索样本空间并相对于 Gibbs 采样减少样本的相关性。为了确保它从后验中采样，它使用先验和似然的乘积来评估提议，这与后验成比例。它根据结果接受的概率来接受，这与当前样本的相应值相关。

提议评估方法的一个关键优点是它使用比后验的确切评估更比例的评估。但是，它可能需要很长时间才能收敛，因为与后验无关的随机移动可能会降低接受率，以至于大量步骤仅产生少量（可能相关的）样本。接受率可以通过减小提议分布的方差来调整，但是由此产生的较小步骤意味着较少的探索。

哈密顿蒙特卡洛 – 走 NUTS

**哈密顿蒙特卡洛（HMC）**是一种混合方法，利用似然梯度的一阶导数信息来提出新的状态以进行探索，并克服了 MCMC 的一些挑战。此外，它还融合了动量以有效地在后验分布中跳跃。因此，与简单的随机游走 Metropolis 或 Gibbs 采样相比，它更快地收敛到高维目标分布。

无转弯采样器是自调节的 HMC 扩展，它自适应地调节在选择提议之前围绕后验的大小和数量。它在高维和复杂的后验分布上表现良好，并且允许拟合许多复杂的模型，而无需对拟合算法本身具有专门的知识。正如我们将在下一节中看到的，它是 PyMC3 中的默认采样器。

变分推断

**变分推断（VI）**是一种通过优化来近似概率密度的机器学习方法。在贝叶斯背景下，它近似后验分布如下：

选择一个参数化的概率分布族
找到该族中距目标最近的成员，以 Kullback-Leibler 散度为度量

与 MCMC 相比，变分贝叶斯往往收敛更快，并且更适用于大型数据。虽然 MCMC 通过链中的样本来近似后验，最终将收敛到任意接近目标，但变分算法通过优化结果来近似后验，这不能保证与目标重合。

变分推断更适合于大型数据集和快速探索多个模型。相比之下，当数据集较小或时间和计算资源的约束较少时，MCMC 将提供更准确的结果。

自动微分变分推断（ADVI）

变分推断的缺点是需要模型特定的导数和实现一个定制的优化例程，这减慢了广泛采用的速度。

最近的**自动微分变分推断（ADVI）**算法自动化了这个过程，用户只需指定模型，以程序形式表达，并且 ADVI 自动生成相应的变分算法（有关实现细节请参阅 GitHub 上的参考资料）。

我们将看到 PyMC3 支持各种变分推断技术，包括 ADVI。

使用 PyMC3 进行概率编程

概率编程提供了一种描述和拟合概率分布的语言，以便我们可以设计、编码和自动估计和评估复杂模型。它旨在抽象掉一些计算和分析复杂性，以使我们能够专注于贝叶斯推理和推断的概念上更为直观和简单的方面。

由于新语言的出现，该领域变得非常动态。Uber 开源了基于 PyTorch 的 Pyro，并且 Google 最近为 TensorFlow 添加了一个概率模块（请参阅 GitHub 上链接的资源）。

结果是，贝叶斯方法在机器学习中的实际相关性和使用可能会增加，以生成关于不确定性的洞见，特别是对于需要透明而不是黑盒模型的用例。

在本节中，我们将介绍流行的 PyMC3 库，该库使用 Python 实现了高级 MCMC 采样和变分推断，用于机器学习模型。与 Stan 一起，以 Monte Carlo 方法的发明者 Stanislaw Ulam 命名，并由哥伦比亚大学的 Andrew Gelman 自 2012 年以来开发，它是最受欢迎的概率编程语言之一。

使用 Theano 的贝叶斯机器学习

PyMC3 于 2017 年 1 月发布，以向 PyMC2（2012 年发布）中使用的 Metropolis-Hastings 采样器添加 Hamiltonian MC 方法。PyMC3 使用 Theano 作为其计算后端，用于动态 C 编译和自动微分。Theano 是一个以矩阵为重点且启用 GPU 的优化库，是由 Yoshua Bengio 的蒙特利尔机器学习算法研究所（MILA）开发的，并受到 TensorFlow 的启发。由于新的深度学习库的成功（有关详细信息，请参阅第十六章《深度学习》），MILA 最近停止进一步开发 Theano。PyMC4 计划于 2019 年使用 TensorFlow，对 API 的影响可能有限。

PyMC3 工作流程

PyMC3 的目标是直观且可读，但功能强大的语法，反映了统计学家描述模型的方式。建模过程通常遵循以下五个步骤：

通过定义以下内容来编码概率模型：
1. 量化关于潜变量的知识和不确定性的先验分布
2. 条件参数在观察数据上的似然函数
使用上一节中描述的选项之一分析后验：
1. 使用 MAP 推断获得一个点估计
2. 使用 MCMC 方法从后验中采样
使用变分贝叶斯近似后验。
使用各种诊断工具检查您的模型。
生成预测。

生成的模型可用于推断，以获取参数值的详细洞察，以及预测新数据点的结果。

我们将使用简单的 logistic 回归来说明这个工作流程（参见笔记本 bayesian_logistic_regression）。随后，我们将使用 PyMC3 来计算和比较贝叶斯夏普比率，估计动态配对交易比率，并实现贝叶斯线性时间序列模型。

模型定义 - 贝叶斯 logistic 回归

如第六章所讨论的，机器学习工作流程， logistic 回归估计一组特征与二进制结果之间的线性关系，通过 S 形函数进行调节，以确保模型产生概率。频率方法导致参数的点估计，这些参数测量了每个特征对数据点属于正类的概率的影响，并且置信区间基于关于参数分布的假设。

相比之下，贝叶斯 logistic 回归估计参数本身的后验分布。后验允许对每个参数的贝叶斯可信区间进行更健壮的估计，其优点在于更透明地了解模型的不确定性。

概率程序由观察到的和未观察到的随机变量（RVs）组成。正如我们所讨论的，我们通过似然分布定义观察到的 RVs，通过先验分布定义未观察到的 RVs。PyMC3 包含许多用于此目的的概率分布。

我们将使用一个简单的数据集，使用每年 $50K 的收入门槛将 30,000 个个体分类。此数据集将包含关于年龄、性别、工作小时和教育年限的信息。因此，我们正在使用这些特征对一个人的收入是否超过$ 50K 的概率进行建模。

PyMC3 库使得对 logistic 回归执行近似贝叶斯推断非常简单。 logistic 回归模型根据以下图表左侧的方式对个体 i 基于 k 个特征的高收入的概率进行建模：

我们将使用上下文管理器 with 来定义一个 manual_logistic_model，以便稍后将其作为概率模型参考：

未观察到的截距和两个特征的参数的随机变量使用假设的先验进行表达，该先验假设为正态分布，均值为 0，标准差为 100。
似然函数根据 logistic 回归的规范将参数与数据结合起来。
结果被建模为 Bernoulli 随机变量，其成功概率由似然函数给出：

with pm.Model() as manual_logistic_model:
    # coefficients as rvs with uninformative priors
    intercept = pm.Normal('intercept', 0, sd=100)
    b1 = pm.Normal('beta_1', 0, sd=100)
    b2 = pm.Normal('beta_2', 0, sd=100)

    # Likelihood transforms rvs into probabilities p(y=1)
    # according to logistic regression model.    
    likelihood = pm.invlogit(intercept + b1 * data.hours + b2 * data.educ)

    # Outcome as Bernoulli rv with success probability 
    # given by sigmoid function conditioned on actual data 
    pm.Bernoulli(name='logit', p=likelihood, observed=data.income)

可视化和板符号

命令 pm.model_to_graphviz(manual_logistic_model) 生成在右侧图中显示的 plate 符号。它显示未观察到的参数为浅色，观察到的元素为深色圆圈。矩形表示由模型定义中包含的数据暗示的观察模型元素的重复次数。

广义线性模型模块

PyMC3 包含许多常见的模型，因此我们通常可以留下自定义应用程序的手动规范。以下代码使用受统计语言 R 启发的公式格式，并由 patsy 库移植到 Python，将相同的逻辑回归定义为 广义线性模型 (GLM) 家族的成员：

with pm.Model() as logistic_model:  
    pm.glm.GLM.from_formula('income ~ hours + educ', 
                            data, 
                            family=pm.glm.families.Binomial())

MAP 推断

我们使用刚刚定义的模型的 .find_MAP() 方法为三个参数获得点 MAP 估计值：

with logistic_model:
    map_estimate = pm.find_MAP()
print_map(map_estimate)
Intercept   -6.561862
hours        0.040681
educ         0.350390

PyMC3 使用拟牛顿 Broyden-Fletcher-Goldfarb-Shanno (BFGS) 算法解决了找到具有最高密度的后验点的优化问题，但提供了几种替代方案，这些替代方案由 sciPy 库提供。结果几乎与相应的 statsmodels 估计相同（有关更多信息，请参阅笔记本）。

近似推断 - MCMC

我们将使用稍微复杂的模型来说明马尔可夫链蒙特卡洛推断：

formula = 'income ~ sex + age+ I(age ** 2) + hours + educ'

Patsy 的函数 I() 允许我们使用常规 Python 表达式动态创建新变量。在这里，我们将 age 平方以捕获更多经验在生活后期增加收入的非线性关系。

请注意，测量尺度非常不同的变量可能会减慢采样过程。因此，我们首先对 age、hours 和 educ 变量进行标准化，应用 sklearn 的 scale() 函数。

一旦我们用新公式定义了我们的模型，我们就可以执行推断以近似后验分布。通过 pm.sample() 函数可用 MCMC 采样算法。

默认情况下，PyMC3 自动选择最有效的采样器，并初始化采样过程以实现有效的收敛。对于连续模型，PyMC3 选择我们在前一节中讨论的 NUTS 采样器。它还通过 ADVI 运行变分推断，以找到采样器的良好起始参数。其中一个选择是使用 MAP 估计。

为了查看收敛情况，我们首先在调整了采样器 1000 次迭代后仅绘制 100 个样本。此后将丢弃这些样本。采样过程可以使用 cores 参数并行化多个链（除非使用 GPU）：

with logistic_model:
    trace = pm.sample(draws=100, tune=1000,
                      init='adapt_diag', # alternative initialization
                      chains=4, cores=2,
                      random_seed=42)

结果跟踪包含每个随机变量的采样值。我们可以通过提供先前运行的跟踪作为输入来继续采样（有关更多信息，请参阅笔记本）。

置信区间

我们可以计算可信区间—贝叶斯的置信区间—作为跟踪的百分位数。结果的边界反映了对于给定概率阈值的参数值范围的信心，而不是参数将在大量试验中多少次在此范围内的数量。笔记本演示了计算和可视化。

近似推断 - 变分贝叶斯

变分推断的界面与 MCMC 实现非常相似。我们只需使用fit()函数而不是sample()函数，还可以选择包括一个早期停止的CheckParametersConvergence回调，如果分布拟合过程收敛到给定的容差：

with logistic_model:
    callback = CheckParametersConvergence(diff='absolute')
    approx = pm.fit(n=100000, 
                    callbacks=[callback])

我们可以从近似分布中抽取样本，以获得类似于之前对 MCMC 采样器所做的跟踪对象：

trace_advi = approx.sample(10000)

检查跟踪摘要显示结果略微不准确。

模型诊断

贝叶斯模型诊断包括验证采样过程是否收敛，并且始终从后验的高概率区域中采样，并确认模型是否很好地代表了数据。

收敛

我们可以随时间和它们的分布可视化样本，以检查结果的质量。下面的图表显示了初始 100 和额外的 100,000 个样本后的后验分布，并说明了收敛意味着多个链识别相同的分布。pm.trace_plot()函数也显示了样本的演变（更多信息请参见笔记本）：

后验分布

PyMC3 为采样器生成各种摘要统计信息。这些信息可以作为 stats 模块中的各个函数提供，或者通过将跟踪提供给pm.summary()函数获取：

	statsmodels	均值	标准差	hpd_2.5	hpd_97.5	n_eff	Rhat
截距	-1.97	-1.97	0.04	-2.04	-1.89	69,492.17	1.00
性别[T.男性]	1.20	1.20	0.04	1.12	1.28	72,374.10	1.00
年龄	1.10	1.10	0.03	1.05	1.15	68,446.73	1.00
I(年龄 ** 2)	-0.54	-0.54	0.02	-0.58	-0.50	66,539.66	1.00
小时	0.32	0.32	0.02	0.28	0.35	93,008.86	1.00
教育	0.84	0.84	0.02	0.80	0.87	98,125.26	1.00

前面的表格中在第一列包括了（分别计算的）statsmodelslogit系数，以显示在这个简单案例中，两个模型是一致的，因为样本均值非常接近系数。

剩下的列包含了最高后验密度（HPD）的估计，用于最小宽度可信区间，这是置信区间的贝叶斯版本，在这里是以 95% 水平计算的。n_eff统计信息总结了 ~100K 绘制结果的有效（未拒绝）样本数量。

R-hat，也称为 Gelman-Rubin 统计量，通过比较链之间的方差与每个链内的方差来检查收敛性。如果采样器收敛，则这些方差应该相同，即链应该看起来相似。因此，统计量应该接近 1。pm.forest_plot() 函数还为多个链总结了此统计量（有关更多信息，请参见笔记本）。

对于具有许多变量的高维模型，检查大量轨迹变得繁琐。使用 NUTS 时，能量图有助于评估收敛问题。它总结了随机过程如何有效地探索后验分布。图表显示了能量和能量转移矩阵，它们应该是相匹配的，如下面的示例所示（有关概念细节，请参见参考资料）：

后验预测检查

后验预测检查（PPCs）非常有用，用于检查模型与数据的拟合程度。它们通过使用来自后验分布的参数生成模型数据来实现此目的。我们使用 pm.sample_ppc 函数进行此操作，并为每个观测值获取 n 个样本（GLM 模块自动将结果命名为 'y'）：

ppc = pm.sample_ppc(trace_NUTS, samples=500, model=logistic_model)
ppc['y'].shape
(500, 29170)

我们可以使用 auc 分数来评估样本内拟合，例如，比较不同模型：

roc_auc_score(y_score=np.mean(ppc['y'], axis=0), 
              y_true=data.income)
0.8294958565103577

预测

在运行后验预测检查之前，预测使用 Theano 的共享变量将训练数据替换为测试数据。为了方便可视化，我们创建一个带有单个预测器小时的变量，创建训练和测试数据集，并将前者转换为共享变量。请注意，我们需要使用 numPy 数组并提供列标签列表（有关详细信息，请参见笔记本）：

X_shared = theano.shared(X_train.values
with pm.Model() as logistic_model_pred:
    pm.glm.GLM(x=X_shared, labels=labels,
               y=y_train, family=pm.glm.families.Binomial())

然后我们像之前一样运行采样器，并在用测试数据替换训练数据后对结果的迹线应用 pm.sample_ppc 函数：

X_shared.set_value(X_test)
ppc = pm.sample_ppc(pred_trace, model=logistic_model_pred,
                    samples=100)

单特征模型的 AUC 分数为 0.65。以下图表显示了每个采样预测器值的实际结果和预测周围的不确定性：

我们现在将说明如何将贝叶斯分析应用于与交易相关的用例。

实际应用

贝叶斯机器学习方法在投资领域有许多应用。概率估计产生的透明度对风险管理和绩效评估自然非常有用。我们将说明如何计算和比较诸如夏普比率之类的指标。GitHub 仓库还包括下面引用的两个笔记本，展示了将贝叶斯 ML 用于建模线性时间序列和随机波动性的用法。

这些笔记本已经改编自 Quantopian 上创建的教程，Thomas Wiecki 领导数据科学，并且在推广贝叶斯方法的使用方面做出了重大贡献。参考资料还包括有关使用贝叶斯 ML 估计配对交易套期保值比率的教程。

贝叶斯夏普比率和表现比较

在本节中，我们将说明如何将夏普比率定义为概率模型，并比较不同收益序列的后验分布。对两组的贝叶斯估计提供了完整的可信值分布，包括效应大小、组均值及其差异、标准差及其差异以及数据的正态性。

主要用例包括分析替代策略之间的差异，或者分析策略的样本内收益与样本外收益之间的差异（详见bayesian_sharpe_ratio笔记本）。贝叶斯夏普比率也是 pyfolio 贝叶斯分析表的一部分。

模型定义

为了将夏普比率建模为概率模型，我们需要关于收益分布和控制此分布的参数的先验。学生 t 分布相对于低自由度（df）的正态分布具有较厚的尾部，是捕捉收益这一方面的合理选择。

因此，我们需要对这个分布的三个参数进行建模，即收益的均值和标准差，以及自由度。我们假设均值和标准差分别服从正态和均匀分布，并且自由度服从具有足够低期望值的指数分布，以确保有厚尾。收益基于这些概率输入，并且年化夏普比率是通过标准计算得出的，忽略了无风险利率（使用每日收益）：

mean_prior = data.stock.mean()
std_prior = data.stock.std()
std_low = std_prior / 1000
std_high = std_prior * 1000

with pm.Model() as sharpe_model:
    mean = pm.Normal('mean', mu=mean_prior, sd=std_prior)
    std = pm.Uniform('std', lower=std_low, upper=std_high)
    nu = pm.Exponential('nu_minus_two', 1 / 29, testval=4) + 2.
    returns = pm.StudentT('returns', nu=nu, mu=mean, sd=std, observed=data.stock)

    sharpe = returns.distribution.mean / returns.distribution.variance ** .5 * np.sqrt(252)
    pm.Deterministic('sharpe', sharpe)

该笔记本包含有关对单个股票进行采样和评估夏普比率的详细信息。

表现比较

为了比较两个收益序列的表现，我们分别对每个组的夏普比率建模，并将效应大小计算为波动率调整后收益之间的差异。通过可视化轨迹，可以深入了解每个指标的分布情况，如下图所示：

用于配对交易的贝叶斯线性回归

在最后一章中，我们介绍了配对交易作为一种流行的算法交易策略，它依赖于两个或更多资产的协整性。给定这样的资产，我们需要估计对冲比率以决定多头和空头仓位的相对大小。基本方法使用线性回归。

linear_regression笔记本说明了贝叶斯线性回归如何跟踪两个资产之间随时间变化的关系。

贝叶斯时间序列模型

PyMC3 包括允许我们对参数不确定性进行类似洞察的 AR(p)模型，与先前的模型相同。bayesian_time_series笔记本说明了一个或多个滞后的时间序列模型。

随机波动模型

正如上一章所讨论的，资产价格具有时变波动性。在某些时期，回报变动很大，而在其他时期则非常稳定。随机波动模型使用潜在波动性变量来建模，该变量被建模为随机过程。无 U 转折采样器是使用这种模型引入的，并且stochastic_volatility笔记本展示了这种用法。

总结

在本章中，我们探讨了机器学习的贝叶斯方法。我们发现它们具有几个优点，包括能够编码先验知识或观点、更深入地了解模型估计和预测周围的不确定性，以及适用于在线学习，在这种情况下，每个训练样本逐渐影响模型的预测。

我们学会了从模型规范到估计、诊断和预测应用贝叶斯工作流程，使用 PyMC3 并探索了几个相关应用。我们将在第十四章中遇到更多贝叶斯模型，主题建模，以及在第十九章中介绍无监督深度学习，在那里我们将介绍变分自动编码器。

接下来的两章介绍基于树的、非线性的集成模型，即随机森林和梯度提升机。

第十章：决策树和随机森林

在本章中，我们将学习两种新的机器学习模型类：决策树和随机森林。我们将看到决策树如何从数据中学习规则，这些规则编码了输入和输出变量之间的非线性关系。我们将说明如何训练决策树并用于回归和分类问题的预测，可视化和解释模型学习到的规则，并调整模型的超参数以优化偏差-方差的权衡并防止过拟合。决策树不仅是重要的独立模型，而且经常被用作其他模型的组成部分。

在本章的第二部分，我们将介绍集成模型，这些模型将多个个体模型组合起来，产生一个具有较低预测误差方差的单一聚合预测。我们将说明自助聚合，通常称为 bagging，作为随机化构建个体模型和减少集成组件预测误差相关性的几种方法之一。

提升是一个非常强大的替代方法，值得拥有自己的章节来讨论一系列最近的发展。我们将说明如何有效地降低方差，并学习如何配置、训练和调整随机森林。我们将看到随机森林作为大量决策树的集合，可以大幅减少预测误差，但会以一定的解释损失为代价。

简而言之，在本章中，我们将涵盖以下内容：

如何使用决策树进行回归和分类
如何从决策树中获得见解，并可视化从数据中学到的决策规则
为什么集成模型往往能够产生更优异的结果
自助聚合如何解决决策树的过拟合挑战
如何训练、调整和解释随机森林

决策树

决策树是一种机器学习算法，它根据从训练数据中学到的决策规则来预测目标变量的值。该算法可以通过改变控制树学习决策规则的目标函数来应用于回归和分类问题。

我们将讨论决策树如何使用规则进行预测，如何训练它们以预测（连续）收益以及（分类）价格走势的方向，以及如何有效地解释、可视化和调整它们。

树如何学习和应用决策规则

我们在第七章和第八章中学习的线性模型通过学习一组参数来预测结果，使用输入变量的线性组合，可能在 logistic 回归的情况下通过 S 形链函数进行转换。

决策树采取了不同的方法：它们学习并顺序地应用一组规则，将数据点分割为子集，然后为每个子集进行一次预测。预测基于应用给定规则序列后产生的训练样本子集的结果值。正如我们将在后面更详细地看到的那样，分类树根据相对类频率或直接的多数类值来预测概率，而回归模型根据可用数据点的结果值均值计算预测。

这些规则中的每一个都依赖于一个特定的特征，并使用阈值将样本分为两组，其值要么低于要么高于该特征的阈值。二叉树自然地表示模型的逻辑：根是所有样本的起点，节点表示决策规则的应用，数据沿着边移动，分割成更小的子集，直到到达一个叶节点，模型在这里做出预测。

对于线性模型，参数值允许解释输入变量对输出和模型预测的影响。相比之下，对于决策树，从根到叶子的路径显示了特征及其值如何通过模型导致特定决策的透明度。

以下图突出显示了模型学习规则的过程。在训练期间，算法扫描特征，并且对于每个特征，试图找到一个分割数据的截断，以最小化由分割产生的损失，加权每个子集中的样本数：

要在训练期间构建整个树，学习算法重复这个划分特征空间的过程，即可能值集合为p个输入变量X[1], X[2], ..., X[p]，划分为互斥且集合完备的区域，每个区域由一个叶节点表示。不幸的是，由于特征和阈值序列的可能组合数量爆炸性增长，算法将无法评估特征空间的每种可能分区。基于树的学习采用自顶向下、贪婪的递归二分拆分方法来克服这种计算限制。

这个过程是递归的，因为它使用了由先前分割产生的数据子集。它是自顶向下的，因为它从树的根节点开始，所有观察结果仍然属于单个区域，然后通过向预测器空间添加一个更多的分割来连续创建树的两个新分支。它是贪婪的，因为算法选择了基于对目标函数的直接影响的最佳规则，而不是展望未来并评估几步之后的损失。我们将在回归和分类树的更具体的上下文中返回到分割逻辑，因为这代表了主要的区别。

随着递归分割向树中添加新节点，训练样本的数量将继续减少。如果规则均匀地分割样本，导致完全平衡的树，每个节点在水平n处将有 2^(n )个节点，每个节点包含总观察数量的相应部分。在实践中，这是不太可能的，因此沿着某些分支的样本数量可能会迅速减少，并且树往往沿着不同路径增长到不同的深度。

为了对新观察结果进行预测，模型使用在训练期间推断出的规则来决定数据点应分配到哪个叶节点，然后使用对应特征空间区域的训练观察结果的均值（用于回归）或众数（用于分类）。特征空间中给定区域的训练样本数较少，即在给定叶节点中，会降低对预测的信心，并可能反映出过度拟合。

递归分割将继续，直到每个叶节点只包含单个样本，并且训练误差已降低到零。我们将介绍几个准则来限制分割并防止决策树产生极端过拟合的自然倾向。

如何在实践中使用决策树

在本节中，我们说明如何使用基于树的模型来获得洞察并进行预测。为了演示回归树，我们预测回报，对于分类情况，我们回到了正面和负面资产价格变动的示例。本节的代码示例在笔记本decision_trees中，除非另有说明。

如何准备数据

我们使用了在 Chapter 4 中构建的数据集的简化版本，Alpha Factor Research。它包含了 Quandl 提供的 2010-2017 年期间的每日股票价格以及各种工程特征。具体详情可以在本章的 GitHub 存储库中的data_prep笔记本中找到。本章中的决策树模型不能处理缺失或分类变量，因此我们将在丢弃任何缺失值后对后者应用虚拟编码。

如何编写自定义交叉验证类

我们还构建了一个针对刚刚创建的数据格式的自定义交叉验证类，该数据具有两个级别的 pandas MultiIndex，一个用于股票代码，另一个用于数据：

class OneStepTimeSeriesSplit:
    """Generates tuples of train_idx, test_idx pairs
    Assumes the index contains a level labeled 'date'"""

    def __init__(self, n_splits=3, test_period_length=1, shuffle=False):
        self.n_splits = n_splits
        self.test_period_length = test_period_length
        self.shuffle = shuffle
        self.test_end = n_splits * test_period_length

    @staticmethod
    def chunks(l, chunk_size):
        for i in range(0, len(l), chunk_size):
            yield l[i:i + chunk_size]

    def split(self, X, y=None, groups=None):
        unique_dates = (X.index
                        .get_level_values('date')
                        .unique()
                        .sort_values(ascending=False)[:self.test_end])

        dates = X.reset_index()[['date']]
        for test_date in self.chunks(unique_dates, self.test_period_length):
            train_idx = dates[dates.date < min(test_date)].index
            test_idx = dates[dates.date.isin(test_date)].index
            if self.shuffle:
                np.random.shuffle(list(train_idx))
            yield train_idx, test_idx

OneStepTimeSeriesSplit 确保了训练和验证集的分割，通过仅使用每只股票的数据直到T-1期进行训练，当使用T月份的数据进行验证时避免了前瞻性偏差。我们将只使用一步预测。

如何构建回归树

回归树基于分配给给定节点的训练样本的平均结果值进行预测，并且通常依赖于均方误差来在递归二分割过程中选择最佳规则。

给定一个训练集，该算法迭代遍历预测变量 X[1], X[2], ..., X[p] 和可能的切分点 s[1], s[1], ..., s[N]，以找到最佳组合。最佳规则将特征空间分成两个区域，{X|X[i] < s[j]} 和 {X|X[i] > s[j]**}，其中 X[i] 特征的值要么低于要么高于 s[j] 阈值，以便基于训练子集的预测最大化相对于当前节点的平方残差的减少。

让我们从一个简化的例子开始，以便于可视化，并仅使用两个月的滞后收益来预测以下月份，类似于上一章的 AR(2) 模型：

使用sklearn，配置和训练回归树非常简单：

from sklearn.tree import DecisionTreeRegressor

# configure regression tree
regression_tree = DecisionTreeRegressor(criterion='mse', # default
                                        max_depth=4,     # up to 4 splits
                                        random_state=42)
# Create training data
y = data.returns
X = data.drop('returns', axis=1)
X2 = X.loc[:, ['t-1', 't-2']]

# fit model
regression_tree.fit(X=X2, y=y)

# fit OLS model
ols_model = sm.OLS(endog=y, exog=sm.add_constant(X2)).fit()

OLS 摘要和决策树的前两个层级的可视化揭示了模型之间的显著差异。OLS 模型提供了三个参数用于截距和两个特征，与该模型对 f 函数的线性假设一致。

相反，回归树图显示了前两个层级的每个节点的特征和阈值用于分割数据（请注意，特征可以重复使用），以及当前均方误差（MSE）的值、样本数量和基于这些训练样本的预测值：

回归树图

该树图还突出显示了样本在节点之间的不均匀分布，只经过两次分裂后，样本数量在 28,000 到 49,000 之间变化。

为了进一步说明关于输入变量与输出之间功能形式关系的不同假设，我们可以将当前收益的预测可视化为特征空间的函数，即作为滞后收益值范围的函数。下图显示了线性回归和回归树的当前期收益作为一期和两期前收益的函数：

右侧的线性回归模型结果突显了滞后和当前收益之间的线性关系，而左侧的回归树图表说明了特征空间的递归划分中编码的非线性关系。

如何构建分类树

分类树的工作方式与回归版本相同，只是结果的分类性质需要不同的方法来进行预测和测量损失。虽然回归树使用相关训练样本的平均结果来预测分配给叶子节点的观测值的响应，但分类树则使用模式，即相关区域中训练样本中最常见的类别。分类树还可以基于相对类频率生成概率预测。

如何优化节点纯度

在构建分类树时，我们也使用递归二分分裂，但是，我们不是使用减少均方误差来评估决策规则的质量，而是可以使用分类错误率，它简单地是给定（叶子）节点中不属于最常见类别的训练样本的比例。

然而，更倾向于使用替代测量方法，如基尼指数或交叉熵，因为它们对节点纯度的敏感性更高，而不是分类错误率。节点纯度是指节点中单一类别的主导程度。一个只包含属于单一类别结果的样本的节点是纯净的，并且暗示着该特征空间区域的成功分类。对于一个分类结果取K个值，0,1,…,K-1，对于表示特征空间区域的给定节点m，其中p[mk]是节点m中第k类结果的比例，它们的计算如下：

当类别比例接近零或一时，基尼不纯度和交叉熵测量值较小，也就是说，当子节点由于分裂而变得纯净时，它们的值最高，并且在二分类情况下，当类别比例均匀或为 0.5 时，它们的值最高。本节末尾的图表显示了这两个测量值和误分类错误率在比例区间 [0, 1] 内的取值情况。

如何训练分类树

我们现在将使用 80% 的样本进行训练，以预测剩余的 20%，训练、可视化和评估一个具有连续 5 次分割的分类树。我们在这里采取了一种简化说明的快捷方式，并使用内置的 train_test_split，它不会防止前瞻性偏差，作为我们的自定义迭代器。树的配置意味着最多有 2⁵=32 个叶节点，平均情况下在平衡的情况下，将包含超过 4,300 个训练样本。看一下以下代码：

# randomize train-test split
X_train, X_test, y_train, y_test = train_test_split(X, y_binary, test_size=0.2, random_state=42)

# configure & train tree learner
classifier = DecisionTreeClassifier(criterion='gini',
                                    max_depth=5,
                                    random_state=42)
classifier.fit(X=X_train, y=y_train)

# Output:
DecisionTreeClassifier(class_weight=None, criterion='gini', max_depth=5,
            max_features=None, max_leaf_nodes=None,
            min_impurity_decrease=0.0, min_impurity_split=None,
            min_samples_leaf=1, min_samples_split=2,
            min_weight_fraction_leaf=0.0, presort=False, random_state=42,
            splitter='best')

在训练模型后，输出显示了我们将在下一节讨论参数调整时更详细讨论的所有 DecisionTreeClassifier 参数。

如何可视化决策树

您可以使用 graphviz 库来可视化树（请参阅 GitHub 获取安装说明），因为 sklearn 可以输出描述树的 .dot 语言，该语言由该库使用。您可以配置输出以包括特征和类别标签，并限制级别数量以使图表可读，如下所示：

dot_data = export_graphviz(classifier,
                           out_file=None, # opt. save to file and convert to png
                           feature_names=X.columns,
                           class_names=['Down', 'Up'],
                           max_depth=3,
                           filled=True,
                           rounded=True,
                           special_characters=True)

graphviz.Source(dot_data)

结果显示，模型使用了各种不同的特征，并指示了连续和分类（虚拟）变量的分割规则。图表显示，在标签值下，每个类别的样本数量，以及在标签类别下，最常见的类别（在样本期间，上升的月份更多）：

如何评估决策树预测结果

为了评估我们的第一个分类树的预测准确性，我们将使用测试集生成预测的类别概率，如下所示：

y_score = classifier.predict_proba(X=X_test)[:, 1] # only keep probabilities for pos. class

.predict_proba() 方法为每个类别生成一个概率。在二分类中，这些概率是互补的，并且总和为 1，因此我们只需要正类的值。为了评估泛化误差，我们将使用基于我们在第六章介绍的接收器操作特性的曲线下面积。结果表明，与基线值 0.5（随机预测）相比，存在显著的改进：

roc_auc_score(y_score=y_score, y_true=y_test)
0.5941

特征重要性

决策树不仅可以被可视化以检查给定特征的决策路径，还可以提供每个特征对拟合到训练数据的模型的贡献的摘要度量。

特征重要性捕获了特征产生的分裂有助于优化模型用于评估分裂质量的度量标准，我们的情况下是基尼不纯度指数。特征的重要性计算为该度量的（归一化）总减少量，并考虑到受分裂影响的样本数量。因此，在树的较早节点使用的特征，其中节点倾向于包含更多样本，通常被认为具有更高的重要性。

下图显示了前 15 个特征的重要性：

过拟合和正则化

决策树有过拟合的倾向，特别是当数据集相对于样本数量具有大量特征时。正如前几章所讨论的，过拟合会增加预测误差，因为模型不仅学习了训练数据中包含的信号，还学习了噪声。

有几种方法可以解决过拟合的风险：

降维（第十二章，无监督学习）通过用更少、更具信息性和更少噪声的特征代表现有特征来改善特征与样本的比例。
集成模型，例如随机森林，结合了多个树，同时在树构建过程中进行随机化，这将在本章的第二部分中介绍。
决策树提供了几个正则化超参数来限制树的增长和相关的复杂性。每个分裂增加节点数，但也减少了每个节点可用于支持预测的样本数量。对于每个额外级别，需要两倍的样本来填充具有相同样本密度的新节点。
树修剪是减少树复杂性的附加工具，通过消除添加了很少价值但增加了模型方差的节点或整个树的部分来实现。例如，成本复杂度剪枝从一个大树开始，通过用叶子替换节点来递归地减小其大小，基本上是将树构建过程反向运行。各个步骤产生的树序列然后可以使用交叉验证进行比较，以选择理想大小。

如何正则化决策树

以下表格列出了 sklearn 决策树实现中可用于此目的的关键参数。在介绍了最重要的参数之后，我们将说明如何使用交叉验证来优化超参数设置，以达到偏差-方差的折衷和降低预测误差：

<tdDefault

参数	选项	描述
`max_depth`	None	int	最大级别数：分割节点直到达到`max_depth`或所有叶子都是纯的或包含少于`min_samples_split`个样本为止。

| **max_features** | None | None：所有特征；int float：分数

auto, sqrt: sqrt(n_features)

log2: log2(n_features) | 用于分割的特征数量。 |

`max_leaf_nodes`	None	None：叶节点数量无限制 int	分割节点直到创建这么多个叶子。
`min_impurity_decrease`	0	float	如果不纯度减少至少这个值，则分割节点。
`min_samples_leaf`	1	int;float (as a percentage of N)	必须在叶节点处的最小样本数。仅当每个左右分支中至少有 `min_samples_leaf` 训练样本时才考虑分割。对于回归问题，可能会平滑模型。
`min_samples_split`	2	int; float (percent of N)	分割内部节点所需的最小样本数：
`min_weight_fraction_leaf`	0		叶节点所需的所有样本权重总和的最小加权分数。除非在 fit 方法中提供了 `sample_weight`，否则样本权重相等。

max_depth 参数对连续分割数量施加硬性限制，并且代表了限制树生长的最直接方法。

min_samples_split 和 min_samples_leaf 参数是限制树生长的替代、数据驱动的方法。与对连续分割数量施加硬性限制不同，这些参数控制进一步分割数据所需的最小样本数。后者保证每个叶子节点有一定数量的样本，而前者在分割导致分布非常不均匀时可能会创建非常小的叶子。小的参数值有助于过度拟合，而高的值可能会阻止树学习数据中的信号。默认值通常相当低，您应该使用交叉验证来探索一系列潜在值。您还可以使用浮点数表示百分比，而不是绝对数值。

sklearn 文档包含有关如何在不同情况下使用各种参数的额外细节；请参阅 GitHub 引用。

决策树剪枝

递归二分法在训练集上可能会产生良好的预测，但往往会过度拟合数据并产生较差的泛化性能，因为它导致了过于复杂的树，表现为大量的叶节点或特征空间的分割。较少的分割和叶节点意味着总体上更小的树，通常也会导致更好的预测性能以及可解释性。

限制叶节点数量的一种方法是避免进一步分割，除非它们对目标度量的改进显著。然而，这种策略的缺点是，有时产生微小改进的分割会在样本的组成不断变化时为后续更有价值的分割创造条件。

相比之下，树剪枝首先通过生长一个非常大的树，然后移除或修剪节点以将大树减少到一个不太复杂且过度拟合的子树。成本复杂度剪枝通过对添加叶节点到树模型的惩罚和调节影响惩罚的正则化参数（类似于套索和岭线性回归模型）生成一系列子树。应用于大树，增加的惩罚将自动产生一系列子树。通过正则化参数的交叉验证可以用来识别最优的修剪子树。

这种方法在 sklearn 中尚不可用；有关进一步详细信息和手动实现剪枝的方法，请参见 GitHub 上的参考资料。

如何调整超参数

决策树提供了一系列超参数来控制和调整训练结果。交叉验证是获取泛化误差无偏估计的最重要工具，这反过来又允许在各种配置选项之间做出知情选择。sklearn 提供了几个工具来简化交叉验证众多参数设置的过程，即我们将在下一节中介绍的GridSearchCV便利类。学习曲线还允许进行诊断，评估收集额外数据以减少泛化误差的潜在益处。

决策树的网格搜索 CV

sklearn提供了一种方法来定义多个超参数的值范围。它自动化了交叉验证这些参数值的各种组合以确定最佳配置的过程。让我们走一遍自动调整模型的过程。

第一步是实例化一个模型对象，并定义一个字典，其中关键词命名超参数，值列出要测试的参数设置：

clf = DecisionTreeClassifier(random_state=42)
param_grid = {'max_depth': range(10, 20),
              'min_samples_leaf': [250, 500, 750],
              'max_features': ['sqrt', 'auto']
              }

然后，实例化GridSearchCV对象，提供评估器对象和参数网格，以及评分方法和交叉验证选择给初始化方法。我们将使用我们自定义的OneStepTimeSeriesSplit类的对象，初始化为使用十个折叠的cv参数，并将评分设置为roc_auc度量。我们可以使用n_jobs参数并行搜索，并通过设置refit=True自动获得使用最佳超参数的训练模型。

所有设置就绪后，我们可以像任何其他模型一样拟合GridSearchCV：

gridsearch_clf = GridSearchCV(estimator=clf,
                          param_grid=param_grid,
                          scoring='roc_auc',
                          n_jobs=-1,
                          cv=cv,  # custom OneStepTimeSeriesSplit
                          refit=True,
                          return_train_score=True)

gridsearch_clf.fit(X=X, y=y_binary)

训练过程为我们的GridSearchCV对象生成一些新属性，最重要的是关于最佳设置和最佳交叉验证分数的信息（现在使用正确的设置以避免前瞻性偏差）。

将max_depth设置为13，min_samples_leaf设置为500，并且在决定分裂时仅随机选择与总特征数的平方根相对应的数量，可以产生最佳结果，AUC 为0.5855：

gridsearch_clf.best_params_
{'max_depth': 13, 'max_features': 'sqrt', 'min_samples_leaf': 500}

gridsearch_clf.best_score_
0.5855

自动化非常方便，但我们也想检查性能如何随不同参数值的变化而变化。完成此过程后，GridSearchCV对象会提供详细的交叉验证结果，以获得更多见解。

如何检查树结构

笔记本还说明了如何手动运行交叉验证以获得自定义树属性，例如与某些超参数设置相关联的总节点数或叶节点数。以下函数访问内部的.tree_属性，以检索有关总节点数以及其中多少节点是叶节点的信息：

def get_leaves_count(tree):
    t = tree.tree_
    n = t.node_count
    leaves = len([i for i in range(t.node_count) if t.children_left[i]== -1])
    return leaves

我们可以将这些信息与训练和测试分数结合起来，以获得有关模型在整个交叉验证过程中行为的详细知识，如下所示：

train_scores, val_scores, leaves = {}, {}, {}
for max_depth in range(1, 26):
    print(max_depth, end=' ', flush=True)
    clf = DecisionTreeClassifier(criterion='gini', 
                                 max_depth=max_depth,
                                 min_samples_leaf=500,
                                 max_features='auto',
                                 random_state=42)
    train_scores[max_depth], val_scores[max_depth], leaves[max_depth] = [], [], []
    for train_idx, test_idx in cv.split(X):
        X_train, y_train,  = X.iloc[train_idx], y_binary.iloc[train_idx]
        X_test, y_test = X.iloc[test_idx], y_binary.iloc[test_idx]
        clf.fit(X=X_train, y=y_train)

        train_pred = clf.predict_proba(X=X_train)[:, 1]
        train_score = roc_auc_score(y_score=train_pred, y_true=y_train)
        train_scores[max_depth].append(train_score)

        test_pred = clf.predict_proba(X=X_test)[:, 1]
        val_score = roc_auc_score(y_score=test_pred, y_true=y_test)
        val_scores[max_depth].append(val_score)    
        leaves[max_depth].append(get_leaves_count(clf))

结果显示在以下图表的左侧面板上。它突出显示了在max_depth设置范围内的样本内外性能，以及围绕误差指标的置信区间。它还显示了在 13 次连续分割中最佳表现的设置，如垂直黑线所示。

学习曲线

学习曲线是一种有用的工具，显示验证和训练分数随训练样本数量的变化而演变。

学习曲线的目的是找出模型是否以及在多大程度上会因为在训练过程中使用更多数据而受益。它还有助于诊断模型的泛化误差更可能是由偏差还是方差驱动。

例如，如果验证分数和训练分数都趋于类似低的值，尽管训练集大小增加，但错误更可能是由于偏差，而额外的训练数据不太可能有所帮助。

请看以下可视化结果：

决策树的优势和劣势

当与迄今为止我们所探索的线性模型相比，回归和分类树在预测时采用非常不同的方法。您如何确定哪个模型更适合手头的问题？请考虑以下内容：

如果结果和特征之间的关系大致是线性的（或者可以相应地进行转换），那么线性回归可能会优于更复杂的方法，例如不利用这种线性结构的决策树。
如果关系呈现高度非线性和更复杂，决策树可能会优于经典模型。

决策树具有几个优点，使其非常受欢迎：

它们相当容易理解和解释，部分因为它们可以很容易地可视化，因此更容易让非技术人员理解。决策树也被称为白盒模型，因为它们在如何得出预测方面具有很高的透明度。黑盒模型，如集成和神经网络，可能会提供更好的预测精度，但是决策逻辑往往更难理解和解释。
与对数据做出更强假设或对数据更敏感的模型（如正则化回归）相比，决策树需要更少的数据准备。
一些决策树实现处理分类输入，不需要创建虚拟变量（提高内存效率），并且可以处理缺失值，就像我们将在第十一章中看到的梯度提升机，但这并不适用于 sklearn。
预测速度快，因为它与叶子节点的数量呈对数关系（除非树变得极不平衡）。
可以使用统计测试验证模型并考虑其可靠性（请参阅 GitHub 参考资料）。

决策树也有一些关键的缺点：

决策树内置了对训练集的过度拟合倾向，并产生高泛化误差。解决这个弱点的关键步骤包括修剪（尚未由 sklearn 支持），以及使用前一节中概述的各种早停准则进行正则化。
与决策树相关的是其高方差，这是由于它们能够紧密地适应训练集，因此数据的微小变化可能会导致决策树结构和因此模型生成的预测产生很大波动。解决决策树高方差的关键机制是使用具有低偏差且产生不相关预测误差的随机决策树集成。
决策树学习的贪婪方法基于局部标准进行优化，即减少当前节点的预测误差，并且不能保证全局最优结果。再次，由随机树组成的集成有助于减轻这个问题。
决策树对不平衡的类权重也很敏感，并可能产生偏倚的树。一种选择是对不平衡的类进行过采样或对更频繁的类进行欠采样。通常最好使用类权重并直接调整目标函数。

随机森林

决策树不仅因其透明度和可解释性而有用，而且是更强大的集成模型的基本构建模块，它将许多个体树与策略结合起来，以随机变化其设计，以解决前一节讨论的过拟合和高方差问题。

集成模型

集成学习涉及将多个机器学习模型组合成一个新模型，旨在比任何单个模型做出更好的预测。更具体地说，一个集成将使用一个或多个给定的学习算法训练的多个基础估计器的预测整合起来，以减少这些模型可能单独产生的泛化错误。

要使集成学习实现这个目标，个体模型必须是：

**准确性：**它们胜过一个天真的基准（例如样本均值或类比例）
**独立：**它们的预测是通过不同方式生成的，以产生不同的错误

集成方法是最成功的机器学习算法之一，特别适用于标准的数值数据。大型集成在机器学习竞赛中非常成功，可能由许多不同的个体模型组成，这些模型通过手工或使用另一个机器学习算法组合在一起。

将不同模型的预测结合起来有几个缺点。这些包括降低的可解释性，以及训练、预测和模型维护的更高复杂性和成本。因此，在实践中（除了竞赛之外），从大规模集成中获得的小幅准确性增益可能不值得增加的成本。

根据它们如何优化组成模型并将结果整合为单个集成预测，通常可以区分两组集成方法：

平均方法独立训练几个基础估计器，然后对它们的预测进行平均。如果基础模型没有偏差并且产生的不高度相关的不同预测错误，那么组合预测可能具有更低的方差，并且可能更可靠。这类似于从具有不相关回报的资产构建投资组合，以减少波动性而不牺牲回报。
提升方法相反，是按顺序训练基础估计器，其特定目标是减少组合估计器的偏差。其动机是将几个弱模型组合成一个强大的集合。

在本章剩余部分，我们将专注于自动平均方法，并在第十一章中讨论梯度提升机的提升方法。

如何降低模型方差

我们发现决策树可能会由于方差较高而做出不良预测，这意味着树结构对训练样本的组成非常敏感。我们还看到低方差的模型，例如线性回归，尽管给定特征数量足够的样本不同，但产生的估计相似。

对于给定的一组独立观察结果，每个都具有σ²的方差，样本均值的标准误差由σ/n给出。换句话说，对更多的观察结果进行平均会减少方差。因此，减少模型方差和其泛化误差的自然方法将是从总体中收集许多训练集，对每个数据集训练不同的模型，并对产生的预测进行平均。

在实践中，我们通常没有许多不同的训练集的奢侈条件。这就是装袋，即自助聚合的缩写，发挥作用的地方。装袋是减少机器学习模型方差的通用方法，特别是在应用于决策树时，特别有用和流行。

袋装是指对自助采样进行聚合，自助采样是带替换的随机样本。这样的随机样本与原始数据集具有相同数量的观察结果，但可能由于替换而包含重复项。

装袋提高了预测准确性，但降低了模型的可解释性，因为不再可能可视化树来理解每个特征的重要性。作为一种集成算法，装袋方法对这些自助采样训练给定数量的基本估计器，然后将它们的预测聚合成最终的集合预测。

装袋通过随机化基本估计器的方法，例如，每棵树的生长方式，然后对预测进行平均，从而减少它们的泛化误差，从而降低了基本估计器的方差。它通常是改进给定模型的直接方法，而无需更改底层算法。它在具有低偏差和高方差的复杂模型中效果最好，例如深度决策树，因为它的目标是限制过拟合。相比之下，提升方法最适合弱模型，例如浅决策树。

有几种装袋方法，它们的不同之处在于它们对训练集应用的随机抽样过程：

粘贴从训练数据中抽取随机样本而不进行替换，而装袋则进行替换
随机子空间无需替换地从特征（即列）中随机抽样
随机补丁通过随机抽样观察和特征来训练基本估计器

袋装决策树

要将装袋应用于决策树，我们从训练数据中创建自助样本，通过反复采样来训练一个决策树，然后在这些样本中的每一个上创建一个决策树，通过对不同树的预测进行平均来创建一个集合预测。

袋装决策树通常生长较大，即具有许多层和叶子节点，并且不进行修剪，以使每棵树具有低偏差但高方差。然后，平均它们的预测的效果旨在减少它们的方差。研究表明，通过构建组合了数百甚至数千棵树的集合的装袋，可以显着提高预测性能。

为了说明装袋对回归树方差的影响，我们可以使用sklearn提供的BaggingRegressor元估计器。它基于指定抽样策略的参数来训练用户定义的基估计器：

max_samples和max_features控制从行和列中抽取的子集的大小，分别
bootstrap和bootstrap_features确定每个样本是有放回还是无放回抽样

以下示例使用指数函数生成单个DecisionTreeRegressor和包含十棵树的BaggingRegressor集成的训练样本，每棵树都生长十层深。这两个模型都是在随机样本上训练的，并为添加了噪声的实际函数预测结果。

由于我们知道真实函数，我们可以将均方误差分解为偏差、方差和噪声，并根据以下分解比较两个模型的这些分量的相对大小：

对于分别包含 250 和 500 个观测值的 100 个重复随机训练和测试样本，我们发现单个决策树的预测方差几乎是基于自举样本的10个装袋树的预测方差的两倍：

noise = .5  # noise relative to std(y)
noise = y.std() * noise_to_signal

X_test = choice(x, size=test_size, replace=False)

max_depth = 10
n_estimators=10

tree = DecisionTreeRegressor(max_depth=max_depth)
bagged_tree = BaggingRegressor(base_estimator=tree, n_estimators=n_estimators)
learners = {'Decision Tree': tree, 'Bagging Regressor': bagged_tree}

predictions = {k: pd.DataFrame() for k, v in learners.items()}
for i in range(reps):
    X_train = choice(x, train_size)
    y_train = f(X_train) + normal(scale=noise, size=train_size)
    for label, learner in learners.items():
        learner.fit(X=X_train.reshape(-1, 1), y=y_train)
        preds = pd.DataFrame({i: learner.predict(X_test.reshape(-1, 1))}, index=X_test)
        predictions[label] = pd.concat([predictions[label], preds], axis=1)

对于每个模型，以下图表显示了上半部分的平均预测值和平均值周围两个标准差的带状区域，以及下半部分基于真实函数值的偏差-方差-噪声分解：

查看笔记本random_forest以获取实现细节。

如何构建随机森林

随机森林算法通过对由装袋生成的自举样本引入的随机化进行扩展，进一步减少方差并提高预测性能。

除了对每个集成成员使用自举训练数据外，随机森林还会对模型中使用的特征进行随机抽样（不重复）。根据实现方式，随机样本可以针对每棵树或每次分裂进行抽取。因此，算法在学习新规则时面临不同选择，无论是在树的级别上还是每次分裂时。

特征样本的大小对回归和分类树有所不同：

对于分类，样本量通常是特征数量的平方根。
对于回归，可以选择从三分之一到所有特征，并应基于交叉验证进行选择。

以下图示说明了随机森林如何随机化单个树的训练，然后将它们的预测聚合成一个集成预测：

除了训练观察值之外，随机化特征的目标是进一步去相关化个体树的预测误差。并非所有特征都是平等的，少量高度相关的特征将在树构建过程中更频繁和更早地被选择，使得决策树在整个集合中更相似。然而，个体树的泛化误差越不相关，整体方差就会减少。

如何训练和调整随机森林

关键的配置参数包括在如何调整超参数部分介绍的各个决策树的各种超参数。以下表格列出了两个RandomForest类的附加选项：

关键字	默认	描述
`bootstrap`	`True`	训练期间使用自举样本。
`n_estimators`	`10`	森林中的树的数量。
`oob_score`	`False`	使用包外样本来估计在未见数据上的 R²。

bootstrap参数激活了前面提到的装袋算法大纲，这反过来又启用了包外分数(oob_score)的计算，该分数使用未包含在用于训练给定树的自举样本中的样本来估计泛化准确度（有关详细信息，请参见下一节）。

n_estimators参数定义了要作为森林一部分生长的树的数量。更大的森林表现更好，但构建时间也更长。重要的是监视交叉验证错误作为基本学习者数量的函数，以确定预测误差的边际减少何时下降，以及额外训练的成本开始超过收益。

max_features参数控制在学习新的决策规则并分裂节点时可用的随机选择特征子集的大小。较低的值会降低树的相关性，从而降低集成的方差，但可能会增加偏差。对于回归问题，良好的起始值是n_features（训练特征的数量），对于分类问题是sqrt(n_features)，但取决于特征之间的关系，并且应该使用交叉验证进行优化。

随机森林设计为包含深度完全生长的树，可以使用max_depth=None和min_samples_split=2来创建。然而，这些值未必是最优的，特别是对于具有许多样本和因此可能非常深的树的高维数据，这可能会导致非常计算密集和内存密集的情况。

sklearn 提供的 RandomForest 类支持并行训练和预测，通过将 n_jobs 参数设置为要在不同核心上运行的 k 个作业数来实现。值 -1 使用所有可用核心。进程间通信的开销可能限制速度提升的线性，因此 k 个作业可能需要超过单个作业的 1/k 时间。尽管如此，在数据庞大且分割评估变得昂贵时，对于大型森林或深度个体树，速度提升通常相当显著，并且可能需要训练相当长时间。

如常，应使用交叉验证来确定最佳参数配置。以下步骤说明了该过程：

我们将使用 GridSearchCV 来识别一组最佳参数，以用于分类树的集成：

rf_clf = RandomForestClassifier(n_estimators=10,
                                criterion='gini',
                                max_depth=None,
                                min_samples_split=2,
                                min_samples_leaf=1,
                                min_weight_fraction_leaf=0.0,
                                max_features='auto',
                                max_leaf_nodes=None,
                                min_impurity_decrease=0.0,
                                min_impurity_split=None,
                                bootstrap=True, oob_score=False,
                                n_jobs=-1, random_state=42)

我们将使用 10 倍自定义交叉验证，并为关键配置设置的值填充参数网格：

cv = OneStepTimeSeriesSplit(n_splits=10)
clf = RandomForestClassifier(random_state=42, n_jobs=-1)
param_grid = {'n_estimators': [200, 400],
              'max_depth': [10, 15, 20],
              'min_samples_leaf': [50, 100]}

使用前述输入配置 GridSearchCV：

gridsearch_clf = GridSearchCV(estimator=clf,
                          param_grid=param_grid,
                          scoring='roc_auc',
                          n_jobs=-1,
                          cv=cv,
                          refit=True,
                          return_train_score=True,
                          verbose=1)

训练由参数网格定义的多个集成模型：

gridsearch_clf.fit(X=X, y=y_binary)

获得最佳参数如下：

gridsearch_clf.bestparams{'max_depth': 15,
 'min_samples_leaf': 100,
 'n_estimators': 400}

最佳分数比单棵树基线略有提高，但具有显著改进：

gridsearch_clf.bestscore_0.6013

随机森林的特征重要性

随机森林集成可能包含数百棵独立的树，但仍然可以从装袋模型中获得关于特征重要性的整体摘要度量。

对于给定的特征，重要性得分是基于该特征进行分割导致的目标函数值的总减少量，平均分布在所有树上。由于目标函数考虑了分割影响的特征数量，所以这个度量隐含地是加权平均的，因此在树的顶部附近使用的特征会由于较少的可用节点中包含的观察次数较多而获得更高的得分。通过对以随机方式生长的许多树进行平均，特征重要性估计失去了一些变化，并且变得更加准确。

基于用于学习决策规则的不同目标，分类和回归树的计算有所不同，分别以回归树的均方误差和分类树的基尼指数或熵来衡量。

sklearn 进一步规范化特征重要性度量，使其总和为 1。因此，计算得出的特征重要性也用作特征选择的一种替代方法，而不是我们在第六章中看到的互信息度量（见 sklearn.feature_selection 模块中的 SelectFromModel）。

在我们的示例中，前 20 个特征的重要性值如下所示：

特征重要性值

包外测试

随机森林提供了内置的交叉验证优势，因为每个树都是在训练数据的自助版本上训练的。因此，每棵树平均使用了仅有的两三分之一可用的观测结果。要了解原因，请考虑自助样本的大小与原始样本的大小相同，每个观测结果被抽取的概率为 1/n。因此，不进入自助样本的概率是(1-1/n)^n，它收敛（迅速）到 1/e，或者大约三分之一。

在用于生长袋装树的训练集中未包括的剩余三分之一的观察结果被称为袋外袋（OOB）观察结果，并且可以作为验证集。正如交叉验证一样，我们对于每棵树都预测未使用此观察结果构建的 OOB 样本的响应，然后对每个 OOB 样本进行单个集合预测的平均预测响应（如果目标是回归）或采取多数投票或预测概率（如果目标是分类）。这些预测产生了关于泛化误差的无偏估计，方便在训练期间计算。

由于预测是在不考虑此观察结果的情况下学习的决策规则产生的，因此结果的 OOB 误差是该观察结果的泛化误差的有效估计。一旦随机森林足够大，OOB 误差就会接近留一交叉验证误差。对于大型数据集，OOB 方法估计测试误差非常高效，而交叉验证可能计算成本高昂。

随机森林的优缺点

袋装集成模型既有优势又有劣势。随机森林的优势包括：

预测性能可以与最佳监督学习算法竞争
它们提供可靠的特征重要性估计
它们提供了测试错误的有效估计，而不需要承担与交叉验证相关的重复模型训练的成本

另一方面，随机森林也有一些缺点：

集合模型本质上比单个决策树不太可解释
训练大量深树可能具有高计算成本（但可以并行化）并且使用大量内存
预测速度较慢，这可能对需要低延迟的应用程序造成挑战

摘要

在本章中，我们学习了一种能够捕捉非线性关系的新模型类，与我们迄今探索过的经典线性模型形成对比。我们看到决策树是如何学习规则来将特征空间划分为产生预测的区域，从而将输入数据分割成特定区域的。

决策树非常有用，因为它们提供了有关特征和目标变量之间关系的独特见解，我们看到了如何可视化树结构中编码的一系列决策规则。

不幸的是，决策树容易过拟合。我们了解到集成模型和自举聚合方法设法克服了决策树的一些缺点，并使它们成为更强大的复合模型的组成部分。

在下一章中，我们将探讨另一个集成模型，它已经被认为是最重要的机器学习算法之一。

第十一章：梯度提升机

在上一章中，我们了解到随机森林通过将它们组合成一个减少个体树高方差的集合来提高个别决策树的预测。随机森林使用装袋（bootstrap aggregation）来将随机元素引入到生长个别树的过程中。

更具体地说，装袋（bagging）从数据中有替换地抽取样本，以便每棵树都在不同但大小相等的数据随机子集上进行训练（其中一些观测重复）。随机森林还随机选择一些特征子集，以便用于训练每棵树的数据的行和列都是原始数据的随机版本。然后，集成通过对个别树的输出进行平均来生成预测。

个别树通常生长较深，以确保低偏差，同时依赖随机化训练过程来产生不同的、不相关的预测错误，这些错误在聚合时具有较低的方差，比个体树的预测更可靠。换句话说，随机化训练旨在使个别树产生的错误彼此不相关或多样化，使集合对过拟合的敏感性大大降低，具有较低的方差，并且对新数据的泛化能力更好。

在本章中，我们将探讨提升（boosting）这种替代机器学习（ML）算法，用于决策树集成，通常能够产生更好的结果。其关键区别在于，提升根据模型在添加新树之前累积的错误来修改用于训练每棵树的数据。与随机森林不同，随机森林独立地使用训练集的不同版本训练许多树，而提升则使用重加权版本的数据进行顺序处理。最先进的提升实现也采用了随机森林的随机化策略。

在本章中，我们将看到提升是如何在过去三十年中演变成最成功的 ML 算法之一的。在撰写本文时，它已经成为结构化数据的机器学习竞赛中的主导者（例如，与高维图像或语音不同，在这些领域中输入和输出之间的关系更加复杂，深度学习表现出色）。具体来说，本章将涵盖以下主题：

提升的工作原理以及与装袋的比较
如何从自适应提升到梯度提升的演变
如何使用 sklearn 使用和调整 AdaBoost 和梯度提升模型
最先进的 GBM 实现如何大幅加速计算
如何防止梯度提升模型过拟合
如何使用xgboost、lightgbm和catboost构建、调整和评估大型数据集上的梯度提升模型
如何解释并从梯度提升模型中获得洞见

自适应提升

像装袋一样，增强是一种集成学习算法，它将基本学习器（通常是决策树）组合成一个集成。增强最初是为分类问题开发的，但也可用于回归，并且已被称为过去 20 年中引入的最有效的学习思想之一（如 Trevor Hastie 等人所述的《统计学习要素》；请查看 GitHub 获取参考链接）。与装袋类似，它是一种通用方法或元方法，可应用于许多统计学习模型。

Boosting 的发展动机是寻找一种方法来将许多弱模型的输出（当一个预测器仅比随机猜测稍好时，称为弱预测器）结合成更强大的，即增强的联合预测。一般来说，Boosting 学习一种形式类似于线性回归的可加假设，H[M]。然而，现在求和的每个元素 m= 1,..., M 都是一个称为 h[t] 的弱基本学习器，它本身需要训练。下面的公式总结了这种方法：

正如上一章所讨论的那样，装袋在不同的训练数据随机样本上训练基本学习器。相比之下，增强是通过在数据上顺序训练基本学习器，该数据反复修改以反映累积学习结果。目标是确保下一个基本学习器补偿当前集合的缺陷。我们将在本章中看到，增强算法在定义缺陷的方式上有所不同。集合使用弱模型的预测的加权平均值进行预测。

第一个带有数学证明，证明它增强了弱学习器性能的提升算法是由罗伯特·舍皮尔（Robert Schapire）和约阿夫·弗洛伊德（Yoav Freund）在 1990 年左右开发的。1997 年，作为分类问题的一个实际解决方案出现了自适应增强（AdaBoost）算法，该算法于 2003 年获得了哥德尔奖。大约又过了五年，该算法在连续梯度下降的方法上得到了推广，利奥·布雷曼（发明了随机森林）将该方法与梯度下降相连接，杰罗姆·弗里德曼在 1999 年提出了梯度增强算法。近年来，出现了许多优化实现，如 XGBoost、LightGBM 和 CatBoost，这些实现牢固地确立了梯度增强作为结构化数据的首选解决方案。

在接下来的几节中，我们将简要介绍 AdaBoost，然后专注于梯度增强模型，以及这个非常强大和灵活的算法的几种最新实现。

AdaBoost 算法

AdaBoost 是第一个在拟合额外集成成员时迭代地适应累积学习进展的增强算法。特别地，AdaBoost 在拟合新的弱学习器之前，会根据当前集成在训练集上的累积错误改变训练数据的权重。AdaBoost 当时是最准确的分类算法，Leo Breiman 在 1996 年的 NIPS 会议上称其为世界上最好的现成分类器。

该算法对 ML 产生了非常重要的影响，因为它提供了理论上的性能保证。这些保证只需要足够的数据和一个可靠地预测略好于随机猜测的弱学习器。由于这种分阶段学习的自适应方法，开发准确的 ML 模型不再需要在整个特征空间上准确地表现。相反，模型的设计可以集中于找到仅仅优于抛硬币的弱学习器。

AdaBoost 与 bagging 有很大的不同，后者构建了在非常深的树上的集成以减少偏差。相比之下，AdaBoost 使用浅树作为弱学习器，通常通过使用树桩来获得更高的准确性，即由单一分裂形成的树。该算法从一个等权重的训练集开始，然后逐步改变样本分布。每次迭代后，AdaBoost 增加被错误分类的观测值的权重，并减少正确预测样本的权重，以便随后的弱学习器更多地关注特别困难的案例。一旦训练完成，新的决策树就会以反映其减少训练错误的贡献的权重被纳入集成中。

对于预测离散类别的一组基本学习器的 AdaBoost 算法，hm，m=1, ..., M，以及N个训练观测值，可以总结如下：

对于观测值 i=1, ..., N，初始化样本权重 w[i]=1/N。
对于每个基本分类器 h[m]，m=1, ..., M，执行以下操作：
1. 用 w[i] 加权训练数据来拟合 hm。
2. 计算训练集上基本学习器的加权错误率 ε[m ]。
3. 计算基本学习器的集成权重 α[m]，作为其错误率的函数，如下式所示：

1. 根据 w[i ] exp(α[m]**)* 更新错误分类样本的权重。
当集成成员的加权和为正时，预测正类；否则，预测负类，如下式所示：

AdaBoost 具有许多实际优势，包括易于实现和快速计算，它可以与任何弱学习器识别方法结合使用。除了集成的大小之外，没有需要调整的超参数。AdaBoost 还对识别异常值很有用，因为接收最高权重的样本是那些一直被错误分类且固有模糊的样本，这也是异常值的典型特征。

另一方面，AdaBoost 在给定数据集上的性能取决于弱学习器充分捕捉特征与结果之间关系的能力。正如理论所示，当数据不足或集成成员的复杂性与数据的复杂性不匹配时，提升将无法表现良好。它也容易受到数据中的噪声影响。

使用 sklearn 的 AdaBoost

作为其集成模块的一部分，sklearn 提供了一个支持两个或多个类的AdaBoostClassifier实现。本节的代码示例位于笔记本gbm_baseline中，该笔记本将各种算法的性能与始终预测最频繁类别的虚拟分类器进行比较。

我们首先需要将base_estimator定义为所有集成成员的模板，然后配置集成本身。我们将使用默认的DecisionTreeClassifier，其中max_depth=1——即一个只有一个分割的树桩。基本估算器的复杂性是关键调参参数，因为它取决于数据的性质。正如前一章所示，对于max_depth的更改应与适当的正则化约束结合使用，例如通过调整min_samples_split来实现，如下代码所示：

base_estimator = DecisionTreeClassifier(criterion='gini', 
                                        splitter='best',
                                        max_depth=1, 
                                        min_samples_split=2, 
                                        min_samples_leaf=20, 
                                        min_weight_fraction_leaf=0.0,
                                        max_features=None, 
                                        random_state=None, 
                                        max_leaf_nodes=None, 
                                        min_impurity_decrease=0.0, 
                                        min_impurity_split=None)

在第二步中，我们将设计集成。n_estimators参数控制弱学习器的数量，learning_rate确定每个弱学习器的贡献，如下代码所示。默认情况下，弱学习器是决策树树桩：

ada_clf = AdaBoostClassifier(base_estimator=base_estimator,
                             n_estimators=200,
                             learning_rate=1.0,
                             algorithm='SAMME.R',
                             random_state=42)

负责良好结果的主要调参参数是n_estimators和基本估算器的复杂性，因为树的深度控制了特征之间的相互作用程度。

我们将使用自定义的 12 折滚动时间序列拆分来交叉验证 AdaBoost 集成，以预测样本中最后 12 个月的 1 个月预测，使用所有可用的先前数据进行训练，如下代码所示：

cv = OneStepTimeSeriesSplit(n_splits=12, test_period_length=1, shuffle=True)
def run_cv(clf, X=X_dummies, y=y, metrics=metrics, cv=cv, fit_params=None):
    return cross_validate(estimator=clf,
                          X=X,
                          y=y,
                          scoring=list(metrics.keys()),
                          cv=cv,
                          return_train_score=True,
                          n_jobs=-1,                    # use all cores
                          verbose=1,
                          fit_params=fit_params)

结果显示加权测试准确度为 0.62，测试 AUC 为 0.6665，负对数损失为-0.6923，以及测试 F1 得分为 0.5876，如下截图所示：

有关交叉验证代码和处理结果的其他详细信息，请参阅配套笔记本。

梯度提升机

AdaBoost 也可以解释为在每次迭代m时通过逐步前向方法最小化指数损失函数来识别与对应权重*α[m]的新基学习器h[m]*添加到集成中，如下式所示：

对 AdaBoost 算法的这种解释是在其发表几年后才发现的。它将 AdaBoost 视为一种基于坐标的梯度下降算法，该算法最小化了特定的损失函数，即指数损失。

梯度提升（Gradient boosting）利用这一洞见，并将提升方法应用于更广泛的损失函数。该方法使得设计机器学习算法以解决任何回归、分类或排名问题成为可能，只要它可以使用可微分的损失函数来表述，并因此具有梯度。定制这一通用方法以适应许多特定预测任务的灵活性对于提升方法的普及至关重要。

梯度提升机（GBM）算法的主要思想是训练基学习器来学习集成当前损失函数的负梯度。因此，每次添加到集成中都直接有助于减少由先前集成成员产生的总体训练误差。由于每个新成员代表数据的新函数，因此也可以说梯度提升是以加法方式优化*h[m]*函数。

简而言之，该算法逐步将弱学习器h[m]（如决策树）拟合到当前集成的损失函数的负梯度上，如下式所示：

换句话说，在给定迭代m时，该算法计算每个观察值的当前损失的梯度，然后将回归树拟合到这些伪残差上。在第二步中，它确定每个终端节点的最佳常数预测，该预测最小化了将该新学习器添加到集成中产生的增量损失。

这与独立决策树和随机森林不同，其中预测取决于相关终端或叶节点中的训练样本的结果值：在回归的情况下是它们的平均值，或者在二元分类的情况下是正类的频率。对损失函数梯度的关注也意味着梯度提升使用回归树来学习回归和分类规则，因为梯度始终是连续函数。

最终集成模型根据个别决策树预测的加权和进行预测，每个决策树都已经训练以在给定一组特征值的情况下最小化集成损失，如下图所示：

梯度提升树在许多分类、回归和排名基准上表现出最先进的性能。它们可能是最受欢迎的集成学习算法，既作为多样化的机器学习竞赛中的独立预测器，也作为实际生产流水线中的一部分，例如，用于预测在线广告的点击率。

梯度提升成功的基础是其以增量方式学习复杂的函数关系的能力。该算法的灵活性需要通过调整约束模型固有倾向于学习训练数据中的噪声而不是信号的超参数来仔细管理过拟合的风险。

我们将介绍控制梯度提升树模型复杂性的关键机制，然后使用 sklearn 实现来说明模型调优。

如何训练和调优 GBM 模型

梯度提升性能的两个关键驱动因素是集成大小和其组成决策树的复杂性。

决策树的复杂性控制旨在避免学习高度具体的规则，这些规则通常意味着叶节点中的样本数量非常少。我们在前一章中介绍了用于限制决策树过拟合到训练数据的最有效约束条件。它们包括要求：

要么分割节点或接受它作为终端节点的最小样本数，或
最小改进节点质量，由纯度或熵或均方误差衡量，对于回归情况而言。

除了直接控制集成大小外，还有各种正则化技术，例如收缩，在第七章 线性模型 中我们遇到的 Ridge 和 Lasso 线性回归模型的上下文中。此外，在随机森林上下文中使用的随机化技术也常用于梯度提升机。

集成大小和提前停止

每次提升迭代旨在减少训练损失，使得对于一个大集成，训练误差可能变得非常小，增加过拟合的风险并在未见数据上表现不佳。交叉验证是找到最小化泛化误差的最佳集成大小的最佳方法，因为它取决于应用程序和可用数据。

由于集成大小需要在训练之前指定，因此监控验证集上的性能并在给定迭代次数时中止训练过程是有用的，当验证误差不再下降时。这种技术称为提前停止，经常用于需要大量迭代并且容易过拟合的模型，包括深度神经网络。

请记住，对于大量试验使用相同的验证集进行早停会导致过拟合，只是过拟合于特定的验证集而不是训练集。最好避免在开发交易策略时运行大量实验，因为误发现的风险显著增加。无论如何，保留一个留置集以获取对泛化误差的无偏估计。

收缩和学习率

收缩技术通过对模型损失函数增加惩罚来应用于增加的模型复杂性。对于提升集合，收缩可以通过将每个新集合成员的贡献缩小一个介于 0 和 1 之间的因子来应用。这个因子被称为提升集合的学习率。降低学习率会增加收缩，因为它降低了每个新决策树对集合的贡献。

学习率与集合大小具有相反的效果，后者往往对较低的学习率增加。较低的学习率与较大的集合结合在一起，已被发现可以降低测试误差，特别是对于回归和概率估计。大量迭代在计算上更昂贵，但通常对于快速的最新实现是可行的，只要个别树保持浅层。根据实现方式，您还可以使用自适应学习率，它会根据迭代次数调整，通常降低后期添加的树的影响。我们将在本章后面看到一些示例。

子抽样和随机梯度提升

如前一章节详细讨论的，自助平均法（bagging）提高了原本噪声分类器的性能。

随机梯度提升在每次迭代中使用无替换抽样，以在训练样本的子集上生成下一棵树。好处是降低计算量，通常可以获得更好的准确性，但子抽样应与收缩结合使用。

正如您所见，超参数的数量不断增加，从而增加了在有限的训练数据上从大量参数试验中选择最佳模型时出现假阳性的风险。最好的方法是依次进行并逐个选择参数值，或者使用低基数子集的组合。

如何在 sklearn 中使用梯度提升

sklearn 的集合模块包含了梯度提升树的实现，用于回归和分类，二元和多类。以下 GradientBoostingClassifier 初始化代码展示了我们之前介绍的关键调整参数，除了我们从独立决策树模型中了解的那些外。笔记本 gbm_tuning_with_sklearn 包含了本节的代码示例。

可用的损失函数包括导致 AdaBoost 算法的指数损失和对应于概率输出的 logistic 回归的 deviance。friedman_mse 节点质量度量是平均平方误差的变化，其中包括一个改进分数（请参阅 GitHub 引用以获取原始论文链接），如以下代码所示：

gb_clf = GradientBoostingClassifier(loss='deviance',                # deviance = logistic reg; exponential: AdaBoost
                                    learning_rate=0.1,              # shrinks the contribution of each tree
                                    n_estimators=100,               # number of boosting stages
                                    subsample=1.0,                  # fraction of samples used t fit base learners
                                    criterion='friedman_mse',       # measures the quality of a split
                                    min_samples_split=2,            
                                    min_samples_leaf=1, 
                                    min_weight_fraction_leaf=0.0,   # min. fraction of sum of weights
                                    max_depth=3,                    # opt value depends on interaction
                                    min_impurity_decrease=0.0, 
                                    min_impurity_split=None, 
                                    max_features=None, 
                                    max_leaf_nodes=None, 
                                    warm_start=False, 
                                    presort='auto',
                                    validation_fraction=0.1, 
                                    tol=0.0001)

与 AdaBoostClassifier 类似，该模型无法处理缺失值。我们将再次使用 12 折交叉验证来获取用于分类滚动 1 个月持有期的方向回报的错误，如下图所示：

gb_cv_result = run_cv(gb_clf, y=y_clean, X=X_dummies_clean)
gb_result = stack_results(gb_cv_result)

我们将解析和绘制结果，以发现与 AdaBoostClassifier 相比略有改进的情况，如下图所示：

如何使用 `GridSearchCV` 调整参数

model_selection 模块中的 GridSearchCV 类促进了对我们想要测试的所有超参数值组合的系统评估。在下面的代码中，我们将说明此功能，用于七个调整参数，当定义时将产生总共 2⁴ x 3² x 4 = 576 种不同的模型配置：

cv = OneStepTimeSeriesSplit(n_splits=12)

param_grid = dict(
        n_estimators=[100, 300],
        learning_rate=[.01, .1, .2],
        max_depth=list(range(3, 13, 3)),
        subsample=[.8, 1],
        min_samples_split=[10, 50],
        min_impurity_decrease=[0, .01],
        max_features=['sqrt', .8, 1]
)

.fit() 方法使用自定义的 OneStepTimeSeriesSplit 和 roc_auc 分数执行交叉验证以评估 12 折。Sklearn 允许我们使用 joblib pickle 实现持久化结果，就像对任何其他模型一样，如以下代码所示：

gs = GridSearchCV(gb_clf,
                   param_grid,
                   cv=cv,
                   scoring='roc_auc',
                   verbose=3,
                   n_jobs=-1,
                   return_train_score=True)
gs.fit(X=X, y=y)

# persist result using joblib for more efficient storage of large numpy arrays
joblib.dump(gs, 'gbm_gridsearch.joblib')

GridSearchCV 对象在完成后具有几个额外的属性，我们可以在加载拾取的结果后访问，以了解哪种超参数组合效果最佳以及其平均交叉验证 AUC 分数，这导致了比默认值略有改进。以下代码显示了这一点：

pd.Series(gridsearch_result.best_params_)
learning_rate              0.01
max_depth                  9.00
max_features               1.00
min_impurity_decrease      0.01
min_samples_split         10.00
n_estimators             300.00
subsample                  0.80

gridsearch_result.best_score_
0.6853

参数对测试分数的影响

GridSearchCV 结果存储了平均交叉验证分数，以便我们分析不同的超参数设置如何影响结果。

下图左侧面板中的六个 seaborn swarm 图显示了所有参数值的 AUC 测试分数分布。在这种情况下，最高的 AUC 测试分数需要较低的 learning_rate 和较大的 max_features 值。某些参数设置，例如较低的 learning_rate，会产生一系列取决于其他参数互补设置的结果范围。其他参数与实验中使用的所有设置的高分数兼容：

我们现在将探讨超参数设置如何共同影响平均交叉验证分数。为了深入了解参数设置是如何相互作用的，我们可以使用DecisionTreeRegressor训练，将平均测试分数作为结果，将参数设置编码为一组独热或虚拟格式的分类变量（详情请参阅笔记本）。树结构突显了使用所有特征（max_features_1），低learning_rate和max_depth超过 3 导致了最佳结果，如下图所示：

在本节第一个图表的右侧面板中的条形图显示了超参数设置对产生不同结果的影响，其特征重要性由决策树产生，决策树生长到其最大深度为止。自然，出现在树顶部附近的特征也累积了最高的重要性得分。

如何在留置集上进行测试

最后，我们希望评估我们从GridSearchCV练习中排除的留置集上最佳模型的性能。它包含了样本期的最后六个月（截至 2018 年 2 月；详情请参阅笔记本）。我们基于 AUC 分数0.6622获得了一个泛化性能估计，使用以下代码：

best_model = gridsearch_result.best_estimator_
preds= best_model.predict(test_feature_data)
roc_auc_score(y_true=test_target, y_score=preds)
0.6622

sklearn 梯度提升实现的缺点是计算速度有限，这使得快速尝试不同的超参数设置变得困难。在接下来的部分中，我们将看到在过去的几年中出现了几种优化实现，大大减少了训练甚至大规模模型所需的时间，并且极大地扩展了这种高效算法的应用范围。

快速可扩展的 GBM 实现

过去几年中，出现了几种新的梯度提升实现，采用了各种创新，加速了训练，提高了资源效率，并允许算法扩展到非常大的数据集。新的实现及其来源如下：

XGBoost（极端梯度提升），2014 年由华盛顿大学的 Tianqi Chen 发起
LightGBM，于 2017 年 1 月由 Microsoft 首次发布
CatBoost，于 2017 年 4 月由 Yandex 首次发布

这些创新解决了训练梯度提升模型的特定挑战（请参阅本章在 GitHub 上的README以获取详细参考）。 XGBoost 的实现是第一个获得流行的新实现：在 2015 年 Kaggle 发布的 29 个获奖解决方案中，有 17 个解决方案使用了 XGBoost。其中 8 个仅依赖于 XGBoost，而其他的则将 XGBoost 与神经网络结合使用。

在说明其实现之前，我们将首先介绍随着时间推移而出现并最终收敛的关键创新（以便大多数功能对所有实现都可用）。

算法创新如何推动性能

随机森林可以通过在独立的自助样本上生长单独的树来并行训练。相比之下，梯度提升的顺序方法会减慢训练速度，这反过来会使得需要适应任务和数据集性质的大量超参数的实验变得更加复杂。

为了通过树扩展整体，训练算法以增量方式最小化相对于整体损失函数的负梯度的预测误差，类似于传统的梯度下降优化器。因此，在训练过程中评估每个特征对决策树与当前梯度拟合的潜在分割点的影响的计算成本与时间成正比。

二阶损失函数近似

最重要的算法创新通过使用依赖于二阶导数的近似来降低评估损失函数的成本，类似于牛顿法找到稳定点。因此，在贪婪树扩展过程中评分潜在分割的速度相对于使用完整损失函数更快。

正如先前提到的，梯度提升模型是以增量方式训练的，其目标是最小化整体 H[M] 的预测误差和正则化惩罚的组合。用 m 步后集成对结果 y[i] 的预测表示为 ŷ[i]^((m))，l 是一个可微的凸损失函数，用来衡量结果和预测之间的差异，Ω 是一个随着整体 H[M] 复杂性增加而增加的惩罚项，增量假设 h[m] 的目标是最小化以下目标函数：

正则化惩罚有助于通过偏好选择使用简单且有预测性的回归树的模型来避免过拟合。例如，在 XGBoost 的情况下，回归树 h 的惩罚取决于每棵树的叶子数 T，每个终端节点的回归树得分 w，以及超参数 γ 和 λ。这在以下公式中总结如下：

因此，在每一步，算法都会贪婪地添加最能改进正则化目标的假设 h[m]。基于泰勒展开的损失函数的二阶近似加速了目标函数的评估，如以下公式所总结：

这里，g[i] 是添加新的学习器前给定特征值的损失函数的一阶梯度，而 h[i] 是相应的二阶梯度（或者海森矩阵）值，如下式所示：

XGBoost 算法是第一个利用损失函数的这种近似来计算给定树结构的最优叶子分数和相应损失函数值的开源算法。分数由终端节点中样本的梯度和 Hessian 的和的比率组成。它使用此值来评分分裂所导致的信息增益，类似于我们在前一章中看到的节点不纯度度量，但适用于任意损失函数（详细推导请参阅 GitHub 上的参考文献）。

简化的分裂查找算法

sklearn 的梯度提升实现找到枚举连续特征所有选项的最优拆分。这种精确的贪婪算法在计算上非常耗费资源，因为它必须首先按特征值对数据进行排序，然后对潜在的非常大数量的拆分选项进行评分和决策。当数据不适合内存或在多台机器上的分布式设置中进行训练时，此方法会面临挑战。

一种近似的分裂查找算法通过将特征值分配给用户确定的一组箱子来减少分裂点的数量，这也可以在训练期间大大减少内存需求，因为每个箱子只需要存储一个分裂。XGBoost 引入了一种分位数草图算法，也能将加权训练样本分成百分位箱子，以实现均匀分布。XGBoost 还引入了处理稀疏数据的能力，这些数据由缺失值、频繁的零梯度统计和一位有效编码引起，并且还可以为给定的分裂学习一个最优的默认方向。因此，该算法只需评估非缺失值。

相比之下，LightGBM 使用基于梯度的单侧采样（GOSS）来排除具有小梯度的大部分样本，并且仅使用其余部分来估算信息增益，并相应地选择分裂值。具有较大梯度的样本需要更多的训练，并且倾向于更多地对信息增益做出贡献。LightGBM 还使用排他性特征绑定来组合那些彼此互斥的特征，它们很少同时取非零值，以减少特征数量。因此，LightGBM 在发布时是最快的实现。

深度优先与叶子优先增长

LightGBM 与 XGBoost 和 CatBoost 在优先拆分哪些节点方面有所不同。LightGBM 以叶子方式决定分裂，即拆分最大化信息增益的叶子节点，即使这会导致不平衡树。相比之下，XGBoost 和 CatBoost 以深度方式扩展所有节点，并在添加更多级别之前首先拆分给定深度的所有节点。这两种方法以不同的顺序扩展节点，除了完整的树之外，它们将产生不同的结果。下图说明了这两种方法：

LightGBM 的叶子优先分割往往增加模型复杂性，并且可能加快收敛速度，但也增加了过拟合的风险。一个深度为 n 的以节点为优先的树最多有 2^(n )个叶子节点，而一个以叶子为优先的树具有 2^n 个叶子，可能会有更多的层次，并且某些叶子节点中包含相对较少的样本。因此，调整 LightGBM 的 num_leaves 设置需要额外的注意，库也允许我们同时控制 max_depth 以避免不必要的节点不平衡。LightGBM 的更近期版本还提供了深度为优先的树增长。

基于 GPU 的训练

所有新的实现支持在一个或多个 GPU 上进行训练和预测，以实现显著的加速。它们与当前的 CUDA 启用的 GPU 兼容。安装要求各不相同，并且正在快速发展。XGBoost 和 CatBoost 的实现适用于几个当前版本，但是 LightGBM 可能需要本地编译（请参阅 GitHub 获取相关文档链接）。

加速取决于库和数据类型，并且范围从低的单位数字倍增到数十倍的因子。仅需更改任务参数并且不需要其他超参数修改就能激活 GPU。

DART – 树的 dropout

2015 年，Rashmi 和 Gilad-Bachrach 提出了一个新模型来训练梯度提升树，旨在解决他们称之为过度专业化的问题：在后续迭代中添加的树往往只影响少数实例的预测，而对于其余实例的贡献较小。然而，该模型的样本外性能可能会受到影响，并且可能对先前在过程中添加的少数树的贡献过于敏感。

新的算法采用了 dropouts，成功地用于学习更准确的深度神经网络，其中 dropouts 在学习过程中静音了部分神经连接。因此，高层节点不能依赖少数连接来传递预测所需的信息。这种方法对于深度神经网络的成功做出了重要贡献，也已与其他学习技术一起使用，例如逻辑回归，以静音特征的随机份额。随机森林和随机梯度提升也会静音随机特征子集。

DART 在树的层次上运作，并且对整棵树进行静音，而不是对个别特征进行操作。使用 DART 生成的集成树的目标是更均匀地对最终预测作出贡献。在某些情况下，这已被证明对排名、回归和分类任务产生更准确的预测。该方法首先在 LightGBM 中实现，并且也适用于 XGBoost。

分类特征的处理

CatBoost 和 LightGBM 实现直接处理分类变量，无需虚拟编码。

CatBoost 的实现（其命名源自对分类特征的处理）包括处理这些特征的几个选项，除了自动独热编码外，还为几个特征的单独分类或多个特征组合分配数字值。换句话说，CatBoost 可以从现有特征的组合中创建新的分类特征。与个别特征或特征组合的分类水平相关的数字值取决于它们与结果值的关系。在分类情况下，这与观察正类别的概率有关，该概率在样本上累积计算，基于先验和平滑系数。有关更详细的数值示例，请参阅文档。

LightGBM 的实现将分类特征的水平分组以最大化（或最小化）与结果值相对于组内的均匀性。

XGBoost 的实现不直接处理分类特征，需要进行独热（或虚拟）编码。

其他功能和优化

XGBoost 通过在内存中保留压缩的列块来优化计算，在几个方面使计算多线程化，其中每个列都按相应特征值排序。XGBoost 在训练前仅计算一次此输入数据布局，并在整个过程中重复使用它以摊销额外的前期成本。使用可以并行完成的分位数时，对列的拆分统计信息的搜索变为线性扫描，并且易于支持列子抽样。

后来发布的 LightGBM 和 CatBoost 库构建在这些创新之上，并通过优化线程和减少内存使用量进一步加速了训练。由于它们的开源性质，这些库随着时间的推移往往会趋于一致。

XGBoost 还支持单调性约束。这些约束确保给定特征的值仅在其整个范围内与结果呈正相关或负相关。它们对于合并已知为真的模型的外部假设非常有用。

如何使用 XGBoost、LightGBM 和 CatBoost

XGBoost、LightGBM 和 CatBoost 提供多种语言的接口，包括 Python，并且具有与其他 sklearn 特性兼容的 sklearn 接口，例如 GridSearchCV，以及它们自己的方法来训练和预测梯度提升模型。gbm_baseline.ipynb 笔记本展示了每个实现的 sklearn 接口的使用。这些库方法通常文档更好，而且也更容易使用，因此我们将使用它们来说明这些模型的使用。

该过程涉及创建特定于库的数据格式，调整各种超参数以及评估我们将在接下来的章节中描述的结果。附带的笔记本包含gbm_tuning.py，gbm_utils.py和gbm_params.py文件，共同提供以下功能，并生成相应的结果。

如何创建二进制数据格式

所有库都有自己的数据格式，用于预先计算特征统计信息以加速搜索分割点，如前所述。这些也可以持久化以加速后续训练的开始。

以下代码为要与OneStepTimeSeriesSplit一起使用的每个模型构造了二进制训练和验证数据集：

cat_cols = ['year', 'month', 'age', 'msize', 'sector']
data = {}
for fold, (train_idx, test_idx) in enumerate(kfold.split(features)):
    print(fold, end=' ', flush=True)
    if model == 'xgboost':
        data[fold] = {'train': xgb.DMatrix(label=target.iloc[train_idx],
                                           data=features.iloc[train_idx],
                                           nthread=-1),                  # use avail. threads
                      'valid': xgb.DMatrix(label=target.iloc[test_idx],
                                           data=features.iloc[test_idx],
                                           nthread=-1)}
    elif model == 'lightgbm':
        train = lgb.Dataset(label=target.iloc[train_idx],
                            data=features.iloc[train_idx],
                            categorical_feature=cat_cols,
                            free_raw_data=False)

        # align validation set histograms with training set
        valid = train.create_valid(label=target.iloc[test_idx],
                                   data=features.iloc[test_idx])

        data[fold] = {'train': train.construct(),
                      'valid': valid.construct()}

    elif model == 'catboost':
        # get categorical feature indices
        cat_cols_idx = [features.columns.get_loc(c) for c in cat_cols]
        data[fold] = {'train': Pool(label=target.iloc[train_idx],
                                    data=features.iloc[train_idx],
                                    cat_features=cat_cols_idx),

                      'valid': Pool(label=target.iloc[test_idx],
                                    data=features.iloc[test_idx],
                                    cat_features=cat_cols_idx)}

可用选项略有不同：

xgboost 允许使用所有可用线程
lightgbm明确地将为验证集创建的分位数与训练集对齐
catboost 实现需要使用索引而不是标签来识别特征列

如何调整超参数

许多超参数在gbm_params.py中列出。每个库都有参数设置来：

指定总体目标和学习算法
设计基学习器
应用各种正则化技术
在训练过程中处理提前停止
启用 GPU 或在 CPU 上进行并行化

每个库的文档详细介绍了可能引用相同概念的各种参数，但是在库之间具有不同名称的参数。 GitHub 存储库包含指向突出显示xgboost和lightgbm相应参数的网站的链接。

目标和损失函数

这些库支持几种提升算法，包括树的梯度提升和线性基学习器，以及 LightGBM 和 XGBoost 的 DART。 LightGBM 还支持我们之前描述的 GOSS 算法，以及随机森林。

梯度提升的吸引力在于有效支持任意可微损失函数，每个库都提供了用于回归，分类和排名任务的各种选项。除了选择的损失函数之外，在训练和交叉验证期间还可以使用其他评估指标来监视性能。

学习参数

梯度提升模型通常使用决策树来捕获特征交互，个体树的大小是最重要的调整参数。 XGBoost 和 CatBoost 将max_depth默认设置为 6。相反，LightGBM 使用默认的num_leaves值为 31，这对应于平衡树的五个级别，但对级别的数量没有约束。为了避免过拟合，num_leaves应低于2^(max_depth)。例如，对于表现良好的max_depth值为 7，您应将num_leaves设置为 70–80，而不是 2⁷=128，或者直接限制max_depth。

树的数量或增强迭代次数定义了整体集合的规模。所有库都支持early_stopping，一旦损失函数在给定迭代次数内不再改善，就会终止训练。因此，通常最好设置大量迭代次数，并根据验证集上的预测性能停止训练。

正则化

所有库都实现了对基本学习器的正则化策略，例如样本数量的最小值或分割和叶节点所需的最小信息增益。

它们还支持通过学习速率通过收缩来在整个集合层次上实现正则化，限制新树的贡献。还可以通过回调函数实现自适应学习速率，随着训练的进行而降低学习速率，这在神经网络的上下文中已经成功使用过。此外，梯度提升损失函数可以通过* L1 或 L2 进行正则化，类似于通过修改Ω（ h [m] *）或通过增加添加更多树的惩罚γ来描述的 Ridge 和 Lasso 线性回归模型。

这些库还允许使用装袋或列子抽样来随机化树的生长，用于随机森林，并减少整体方差以去相关预测误差。为了保护免受过拟合，对于近似分割查找的特征量化，添加更大的箱子作为另一种选择。

随机化网格搜索

为了探索超参数空间，我们为我们想要测试的关键参数指定值。sklearn 库支持RandomizedSearchCV，从指定的分布中随机抽样一部分参数组合进行交叉验证。我们将实现一个自定义版本，允许我们利用早停，同时监视当前表现最佳的组合，因此我们可以在满意结果时中止搜索过程，而不是事先指定一组迭代次数。

为此，我们根据每个库的参数指定参数网格，使用itertools库提供的内置笛卡尔product生成器生成所有组合，并随机shuffle结果。在 LightGBM 的情况下，我们会自动根据当前num_leaves值设置max_depth，如下所示：

param_grid = dict(
        # common options
        learning_rate=[.01, .1, .3],
        colsample_bytree=[.8, 1],  # except catboost

        # lightgbm
        num_leaves=[2 ** i for i in range(9, 14)],
        boosting=['gbdt', 'dart'],
        min_gain_to_split=[0, 1, 5],  # not supported on GPU

all_params = list(product(*param_grid.values()))
n_models = len(all_params) # max number of models to cross-validate
shuffle(all_params)

然后，我们执行交叉验证如下：

GBM = 'lightgbm'
for test_param in all_params:
    cv_params = get_params(GBM)
    cv_params.update(dict(zip(param_grid.keys(), test_param)))
    if GBM == 'lightgbm':
        cv_params['max_depth'] = int(ceil(np.log2(cv_params['num_leaves'])))
    results[n] = run_cv(test_params=cv_params,
                        data=datasets,
                        n_splits=n_splits,
                        gb_machine=GBM)

run_cv函数实现了三个库的交叉验证。对于light_gbm示例，该过程如下：

def run_cv(test_params, data, n_splits=10):
    """Train-Validate with early stopping"""
    result = []
    cols = ['rounds', 'train', 'valid']
    for fold in range(n_splits):
        train = data[fold]['train']
        valid = data[fold]['valid']

        scores = {}
        model = lgb.train(params=test_params,
                          train_set=train,
                          valid_sets=[train, valid],
                          valid_names=['train', 'valid'],
                          num_boost_round=250,
                          early_stopping_rounds=25,
                          verbose_eval=50,
                          evals_result=scores)

        result.append([model.current_iteration(),
                       scores['train']['auc'][-1],
                       scores['valid']['auc'][-1]])

    return pd.DataFrame(result, columns=cols)

train()方法还会生成存储在scores字典中的验证分数。当早停生效时，最后一次迭代也是最佳分数。有关更多详细信息，请参见 GitHub 上的完整实现。

如何评估结果

使用 GPU，我们可以在几分钟内训练一个模型，并在几个小时内评估数百个参数组合，而使用 sklearn 实现则需要多天。对于 LightGBM 模型，我们探索了使用库处理分类变量的因子版本和使用独热编码的虚拟版本。

结果存储在 model_tuning.h5 HDF5 存储中。模型评估代码样本在 eval_results.ipynb 笔记本中。

跨模型的交叉验证结果

当比较四次测试运行中三个库的平均交叉验证 AUC 时，我们发现 CatBoost 为表现最佳模型产生了稍高的 AUC 分数，同时也产生了最广泛的结果分布，如下图所示：

表现最佳的 CatBoost 模型使用以下参数（详见笔记本）：

max_depth 为 12，max_bin 为 128
max_ctr_complexity 为 2，限制了分类特征的组合数量
one_hot_max_size 为 2，排除了二元特征的数值变量分配
random_strength 不等于 0 以随机化分裂的评估

训练相对于 LightGBM 和 XGBoost 稍慢（都使用 GPU），平均每个模型 230 秒。

对 LightGBM 和 XGBoost 模型表现最好的更详细的分析显示，LightGBM 因子模型的性能几乎与其他两个模型相当，但模型复杂度要低得多。它平均只包含 41 棵树，深度为三级，每棵树最多有八个叶子节点，并且还使用了 min_gain_to_split 形式的正则化。它在训练集上过拟合明显较少，训练 AUC 仅略高于验证 AUC。它还训练速度更快，每个模型只需 18 秒，因为它的复杂度更低。实际上，这个模型更可取，因为它更有可能产生良好的样本外表现。具体细节如下表所示：

	LightGBM 虚拟	XGBoost 虚拟	LightGBM 因子
验证 AUC	68.57%	68.36%	68.32%
训练 AUC	82.35%	79.81%	72.12%
`learning_rate`	0.1	0.1	0.3
`max_depth`	13	9	3
`num_leaves`	8192		8
`colsample_bytree`	0.8	1	1
`min_gain_to_split`	0	1	0
轮数	44.42	59.17	41.00
时间	86.55	85.37	18.78

下图显示了不同 max_depth 设置对 LightGBM 和 XGBoost 模型的验证分数的影响：较浅的树产生了更广泛的结果范围，需要与适当的学习率和正则化设置相结合，以产生前面表格中显示的强结果：

与之前显示的 DecisionTreeRegressor 不同，我们也可以使用线性回归来评估不同特征在验证 AUC 分数方面的统计显著性。对于 LightGBM 虚拟模型，其中回归解释了结果变异的 68%，我们发现只有 min_gain_to_split 正则化参数不显著，如下图所示：

在实践中，深入了解模型如何进行预测非常重要，特别是对于投资策略，决策者通常需要合理的解释。

如何解释 GBM 结果

理解模型为什么预测某个结果对于多种原因都非常重要，包括信任、可操作性、责任和调试。通过模型发现的特征与结果之间的非线性关系以及特征之间的相互作用也在学习更多关于所研究现象的基本驱动因素时具有价值。

了解树集成方法（如梯度提升或随机森林模型）所做预测的见解的常见方法是将特征重要性值归因于每个输入变量。这些特征重要性值可以根据单个预测的情况或整个数据集（即所有样本）全局计算，以获得模型进行预测的更高级别视角。

特征重要性

有三种主要方法来计算全局特征重要性值：

Gain: 这种经典方法由 Leo Breiman 在 1984 年引入，使用给定特征的所有分割所贡献的损失或杂质的总减少。其动机在很大程度上是启发式的，但它是一种常用的特征选择方法。
分割计数：这是一种替代方法，它计算使用特征进行分割决策的频率，根据选择特征以达到所得信息增益。
排列：这种方法随机排列测试集中的特征值，并测量模型误差的变化量，假设重要特征应该会导致预测误差大幅增加。不同的排列选择导致该基本方法的替代实现。

对于单个预测计算特征重要性值的方法较少，因为可用的模型无关解释方法比树特定方法慢得多。

所有梯度提升实现在训练后都提供特征重要性得分作为模型属性。XGBoost 库提供五个版本，如下列表所示：

total_gain 和 gain 作为其每个分割的平均值
total_cover 作为使用特征时每个分割的样本数
weight 作为前面值的分割计数

可使用训练模型的.get_score()方法和相应的importance_type参数获取这些值。对于性能最佳的 XGBoost 模型，结果如下（total度量具有 0.8 的相关性，cover和total_cover也是如此）：

虽然不同月份和年份的指标占主导地位，但最近 1 个月的回报是从total_gain的角度来看第二重要的特征，并且根据weight度量经常被使用，但由于它平均应用于相对较少的实例而产生较低的平均收益（有关实施细节，请参见笔记本）。

偏依赖图

除了个体特征对模型预测的总体贡献之外，偏依赖图还可视化目标变量与一组特征之间的关系。梯度提升树的非线性特性导致这种关系取决于所有其他特征的值。因此，我们将对这些特征进行边际化处理。通过这样做，我们可以将偏依赖性解释为期望的目标响应。

我们只能为单个特征或特征对可视化偏依赖性。后者会产生等高线图，显示不同特征值的组合如何产生不同的预测概率，如下所示：

fig, axes = plot_partial_dependence(gbrt=gb_clf,
                                    X=X_dummies_clean,
                                    features=['month_9', 'return_1m', 'return_3m', ('return_1m', 'return_3m')],
                                    feature_names=['month_9','return_1m', 'return_3m'],
                                    percentiles=(0.01, 0.99),
                                    n_jobs=-1,
                                    n_cols=2,
                                    grid_resolution=250)

经过一些额外的格式化（请参阅配套笔记本），我们得到了以下图表：

右下角的图显示了在消除[1%，99%]百分位数处的离群值后，给定滞后 1 个月和 3 个月收益值范围的情况下，下个月产生正回报的概率的依赖性。 month_9 变量是一个虚拟变量，因此呈阶梯函数样式的图。我们还可以将依赖性可视化为 3D，如下代码所示：

targets = ['return_1m', 'return_3m']
target_feature = [X_dummies_clean.columns.get_loc(t) for t in targets]
pdp, axes = partial_dependence(gb_clf,
                               target_feature,
                               X=X_dummies_clean,
                               grid_resolution=100)

XX, YY = np.meshgrid(axes[0], axes[1])
Z = pdp[0].reshape(list(map(np.size, axes))).T

fig = plt.figure(figsize=(14, 8))
ax = Axes3D(fig)
surf = ax.plot_surface(XX, YY, Z,
                       rstride=1,
                       cstride=1,
                       cmap=plt.cm.BuPu,
                       edgecolor='k')
ax.set_xlabel(' '.join(targets[0].split('_')).capitalize())
ax.set_ylabel(' '.join(targets[1].split('_')).capitalize())
ax.set_zlabel('Partial Dependence')
ax.view_init(elev=22, azim=30)

这产生了以下 1 个月回报方向对滞后 1 个月和 3 个月回报的偏依赖的 3D 图：

SHapley Additive exPlanations

在 2017 年 NIPS 会议上，来自华盛顿大学的 Scott Lundberg 和 Su-In Lee 提出了一种新的更准确的方法来解释树集成模型对输出的贡献，称为SHapley Additive exPlanations，或SHAP值。

这个新算法基于这样一个观察：树集成的特征归因方法，如我们之前看过的方法，是不一致的——即，使特征对输出的影响增加的模型变化可能会降低该特征的重要性值（请参阅 GitHub 上有关此的详细说明）。

SHAP 值统一了协作博弈理论和局部解释的思想，并且根据期望值显示其在理论上是最优、一致和局部准确的。最重要的是，Lundberg 和 Lee 开发了一种算法，成功将计算这些与模型无关的、可加的特征归因方法的复杂度从 O(TLD^M) 减少到 O(TLD²)，其中 T 和 M 分别是树和特征的数量，D 和 L 是树的最大深度和叶子节点数。这一重要的创新使得对先前难以处理的具有数千棵树和特征的模型的预测能够在几秒钟内解释。一个开源实现于 2017 年底发布，并且兼容 XGBoost、LightGBM、CatBoost 和 sklearn 树模型。

Shapley 值源自博弈论，是一种为协作游戏中的每个玩家分配值的技术，反映了他们对团队成功的贡献。SHAP 值是对基于树的模型的博弈论概念的一种改编，并且为每个特征和每个样本计算。它们衡量了一个特征对给定观察结果的模型输出的贡献。因此，SHAP 值提供了不同的洞察力，说明了特征的影响如何在样本间变化，这在这些非线性模型中交互作用效应的作用中非常重要。

如何按特征总结 SHAP 值

要获得对多个样本的特征重要性的高级概述，有两种绘制 SHAP 值的方法：对所有样本进行简单平均，类似于先前计算的全局特征重要性度量（如下图左侧面板所示），或者绘制散点图来显示每个样本的每个特征的影响（如下图右侧面板所示）。使用兼容库的训练模型和匹配输入数据非常容易产生，如下面的代码所示：

# load JS visualization code to notebook
shap.initjs()

# explain the model's predictions using SHAP values
explainer = shap.TreeExplainer(model)
shap_values = explainer.shap_values(X_test)

shap.summary_plot(shap_values, X_test, show=False)

下图右侧的散点图按照所有样本的 SHAP 值对特征进行排序，然后显示了每个特征对模型输出的影响，其值由特征的值的函数表示，通过其颜色表示，红色表示高值，蓝色表示低值，相对于特征的范围：

如何使用力图解释预测

下面的力图显示了各种特征及其值对模型输出的累积影响，在这种情况下，模型输出为 0.6，比基准值 0.13（提供的数据集上的平均模型输出）高得多。突出显示为红色的特征增加了输出。十月份是最重要的特征，将输出从 0.338 增加到 0.537，而 2017 年则降低了输出。

因此，我们获得了模型如何到达特定预测的详细分解，如下图所示：

我们还可以同时计算大量数据点或预测的力图，并使用聚类可视化来深入了解数据集中某些影响模式的普遍程度。以下图表显示了前 1000 个观察结果的力图，旋转了 90 度，水平堆叠，并按照不同特征对给定观察结果的影响排序。该实现使用层次凝聚聚类对特征 SHAP 值上的数据点进行标识以识别这些模式，并以交互方式显示结果进行探索性分析（参见笔记本），如以下代码所示：

shap.force_plot(explainer.expected_value, shap_values[:1000,:], X_test.iloc[:1000])

这将产生以下输出：

如何分析特征交互

最后，SHAP 值使我们能够通过将这些相互作用与主要效应分离来获得对不同特征之间的相互作用效应的额外洞察。shap.dependence_plot 可以定义如下：

shap.dependence_plot("return_1m", shap_values, X_test, interaction_index=2, title='Interaction between 1- and 3-Month Returns')

它显示了 1 个月回报的不同值（x 轴）如何影响结果（y 轴上的 SHAP 值），并由 3 个月回报区分：

SHAP 值在每个单独预测的水平上提供了细粒度的特征归因，并通过（交互式）可视化实现了对复杂模型的更丰富检查。本节开头显示的 SHAP 摘要散点图提供了比全局特征重要性条形图更加差异化的见解。单个聚类预测的力图允许进行更详细的分析，而 SHAP 依赖图捕获了相互作用效应，并因此提供了比局部依赖图更准确和详细的结果。

SHAP 值的限制与任何当前特征重要性度量一样，涉及对高度相关的变量影响的归因，因为它们相似的影响可能以任意方式被分解。

摘要

在本章中，我们探讨了梯度提升算法，该算法用于以顺序方式构建集成，添加一个只使用非常少的特征的浅层决策树，以改善已经做出的预测。我们看到了梯度提升树如何可以非常灵活地应用于广泛的损失函数，并提供了许多机会来调整模型以适应给定的数据集和学习任务。

最近的实现大大简化了梯度提升的使用，通过加速训练过程并提供更一致和详细的特征重要性以及单个预测驱动因素的洞察。在下一章中，我们将转向贝叶斯方法来进行机器学习。

Python-机器学习算法交易实用指南-三-