贝叶斯加法回归树：简介#

from pathlib import Path

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import pymc as pm
import pymc_bart as pmb

from sklearn.model_selection import train_test_split

%config InlineBackend.figure_format = "retina"

print(f"Running on PyMC v{pm.__version__}")

Running on PyMC v5.19.1

RANDOM_SEED = 5781
np.random.seed(RANDOM_SEED)
az.style.use("arviz-darkgrid")

BART 概述#

贝叶斯加法回归树 (BART) 是一种非参数回归方法。如果我们有一些协变量 \(X\) 并且我们想用它们来建模 \(Y\)，一个 BART 模型（省略先验）可以表示为

\[Y = f(X) + \epsilon\]

其中我们使用 \(m\) 个 回归树 的总和来建模 \(f\)，并且 \(\epsilon\) 是一些噪声。在最典型的例子中，\(\epsilon\) 是正态分布的，\(\mathcal{N}(0, \sigma)\)。所以我们也可以写成

\[Y \sim \mathcal{N}(\mu=BART(X), \sigma)\]

原则上，没有任何限制我们使用树的总和来建模其他关系。例如，我们可能有

\[Y \sim \text{Poisson}(\mu=BART(X))\]

BART 是贝叶斯方法的原因之一是在回归树上使用了先验。先验的定义方式是它们偏爱叶值接近于零的浅树。一个关键思想是，单个 BART 树在拟合数据方面不是很擅长，但是当我们将许多这些树相加时，我们会得到一个良好且灵活的近似。

使用 BART 进行煤矿开采#

为了更好地理解 BART 在实践中的应用，我们将使用经典的煤矿灾难数据集。这是 PyMC 中的经典示例之一。我们不将此问题视为具有两个泊松分布的切换点模型，就像原始 PyMC 示例中那样。我们将把这个问题看作是具有泊松响应的非参数回归（这通常在 泊松过程 或 Cox 过程 方面进行讨论，但我们可以不用深入研究这些技术细节）。对于类似的示例，但使用高斯过程，请参阅 1 或 2。由于我们的数据只是一个包含日期的单列，我们需要进行一些预处理。我们将对数据进行离散化，就像我们正在构建直方图一样。我们将使用箱子的中心作为变量 \(X\)，并将每个箱子的计数作为变量 \(Y\)

try:
    coal = np.loadtxt(Path("..", "data", "coal.csv"))
except FileNotFoundError:
    coal = np.loadtxt(pm.get_data("coal.csv"))

# discretize data
years = int(coal.max() - coal.min())
bins = years // 4
hist, x_edges = np.histogram(coal, bins=bins)
# compute the location of the centers of the discretized data
x_centers = x_edges[:-1] + (x_edges[1] - x_edges[0]) / 2
# xdata needs to be 2D for BART
x_data = x_centers[:, None]
# express data as the rate number of disaster per year
y_data = hist

在 PyMC 中，BART 变量的定义与其他随机变量非常相似。一个重要的区别是，我们必须将我们的 Xs 和 Ys 传递给 BART 变量，此信息在对树进行采样时使用，树的总和的先验非常大，如果没有来自我们数据的任何信息，这将是一项不可能完成的任务。

在这里，我们还明确指出我们将使用 20 棵树的总和 (m=20)。像 20 这样少的树数量对于像这样的简单模型来说可能已经足够好了，并且作为更复杂模型的快速近似也可能非常有效，尤其是在建模的早期阶段，那时我们可能想要快速尝试一些东西，以便更好地掌握哪个模型可能对我们的问题来说是一个好主意。在这些情况下，一旦我们对我们真正喜欢的模型有更多的确定性，我们可以通过增加 m 来改进近似，在文献中，常见的是找到关于使用 50、100 或 200 等数字获得良好结果的报告。

with pm.Model() as model_coal:
    μ_ = pmb.BART("μ_", X=x_data, Y=np.log(y_data), m=20)
    μ = pm.Deterministic("μ", pm.math.exp(μ_))
    y_pred = pm.Poisson("y_pred", mu=μ, observed=y_data)
    idata_coal = pm.sample(random_seed=RANDOM_SEED)

Multiprocess sampling (4 chains in 4 jobs)
PGBART: [μ_]





Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 11 seconds.
The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details









在检查结果之前，我们需要讨论一个更详细的细节，BART 变量始终在实线上进行采样，这意味着原则上我们可以获得从 \(-\infty\) 到 \(\infty\) 的值。因此，我们可能需要像对标准广义线性模型那样转换它们的值，例如在 model_coal 中，我们计算了 pm.math.exp(μ_)，因为泊松分布期望的值是从 0 到 \(\infty\)。这是通常的做法，新颖之处在于我们可能需要将逆变换应用于 Y 的值，就像我们在之前的模型中所做的那样，我们取了 \(\log(Y)\)。这样做的主要原因是 Y 的值用于获得树总和的合理初始值以及叶节点的方差。因此，应用逆变换是提高结果效率和准确性的一种简单方法。我们是否应该对每种可能的似然函数都这样做？好吧，不。如果我们正在使用 BART 来表示正态分布、学生 t 分布或非对称拉普拉斯分布等分布的位置参数，我们不需要做任何事情，因为这些参数的支持也是实线。一个重要的例外是伯努利似然函数（或 n=1 的二项分布），在这种情况下，我们需要将 logistic 函数应用于 BART 变量，但无需应用其逆变换来转换 Y，PyMC-BART 已经处理了这种情况。


好的，现在让我们看看 model_coal 的结果。






_, ax = plt.subplots(figsize=(10, 6))

rates = idata_coal.posterior["μ"] / 4
rate_mean = rates.mean(dim=["draw", "chain"])
ax.plot(x_centers, rate_mean, "w", lw=3)
ax.plot(x_centers, y_data / 4, "k.")
az.plot_hdi(x_centers, rates, smooth=False)
az.plot_hdi(x_centers, rates, hdi_prob=0.5, smooth=False, plot_kwargs={"alpha": 0})
ax.plot(coal, np.zeros_like(coal) - 0.5, "k|")
ax.set_xlabel("years")
ax.set_ylabel("rate");








../_images/799c9ad4685dcfcd9174007a8904d3fb4151861c2e66f93005878b67d6bdd470.png





下图中的白线显示了事故发生率的中位数。较深的橙色带代表 HDI 50%，较浅的橙色带代表 94%。我们可以看到 1880 年至 1900 年间煤矿事故迅速减少。请随意将这些结果与原始 PyMC 简介概述 示例中的结果进行比较。


在之前的图中，白线是 4000 个后验抽样的平均值，每个后验抽样是 m=20 棵树的总和。


下图显示了 \(\mu\) 的后验的两个样本。我们可以看到这些函数不是平滑的。这很好，并且是使用回归树的直接结果。树可以被看作是表示阶梯函数的一种方式，而阶梯函数的总和只是另一个阶梯函数。因此，当使用 BART 时，我们需要知道我们假设阶梯函数对于我们的问题来说是一个足够好的近似。在实践中，情况通常如此，因为我们对许多树求和，通常值如 50、100 或 200。此外，我们经常对后验分布求平均值。所有这些都使“步进更平滑”，即使我们从来没有真正得到一个平滑函数，例如使用高斯过程或样条。一个很好的理论结果告诉我们，在 \(m \to \infty\) 的极限情况下，BART 先验收敛到 处处不可微 的高斯过程。


下图显示了 \(\mu\) 从后验的两个样本。






plt.step(x_data, rates.sel(chain=0, draw=[3, 10]).T);








../_images/282353251f94d965df18ff88367089a674a42d712285bc39b0eee1570358e28d.png









使用 BART 进行自行车骑行#


为了探索 PyMC-BART 提供的其他功能。我们现在将继续讨论另一个例子。在本例中，我们有关于城市自行车租赁数量的数据，并且我们选择了四个协变量；一天中的 hour、temperature、humidity 以及是否是 workingday 还是周末。此数据集是 bike_sharing_dataset 的子集。






try:
    bikes = pd.read_csv(Path("..", "data", "bikes.csv"))
except FileNotFoundError:
    bikes = pd.read_csv(pm.get_data("bikes.csv"))

features = ["hour", "temperature", "humidity", "workingday"]

X = bikes[features]
Y = bikes["count"]













with pm.Model() as model_bikes:
    α = pm.Exponential("α", 1)
    μ = pmb.BART("μ", X, np.log(Y), m=50)
    y = pm.NegativeBinomial("y", mu=pm.math.exp(μ), alpha=α, observed=Y)
    idata_bikes = pm.sample(compute_convergence_checks=False, random_seed=RANDOM_SEED)









Multiprocess sampling (4 chains in 4 jobs)
CompoundStep
>NUTS: [α]
>PGBART: [μ]









Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 37 seconds.











收敛诊断#


为了检查 BART 模型的采样收敛性，我们建议采用两步法。


    

  • 对于非 BART 变量（如 model_bikes 中的 \(\alpha\)），我们遵循标准建议，例如检查 R-hat (<= 1.01) 和 ESS (< 100x 链数) 数值诊断，以及使用迹线图，甚至更好的秩图
    • 

  • 对于 BART 变量，我们建议使用 pmb.plot_convergence 函数。
    • 



我们可以在下面看到这些检查






az.plot_trace(idata_bikes, var_names=["α"], kind="rank_bars");








../_images/cde16cda654ca41988d8717ab5da2f6dd538f78927ebeb248f77948500d056e9.png









pmb.plot_convergence(idata_bikes, var_name="μ");









/home/osvaldo/proyectos/00_BM/arviz-devs/arviz/arviz/plots/ecdfplot.py:298: BehaviourChangeWarning: In future versions, if `eval_points` is not provided, then the ECDF will be evaluated at the unique values of the sample. To keep the current behavior, provide `eval_points` explicitly.
  warnings.warn(




../_images/45d13c54ac318b12cd4866ea8342fee13a20eebd3d2151d67edfd92bf50a7c7b.png





在 BART 文献中，BART 变量的诊断有时被认为不如非 BART 变量的诊断重要，主要论点是潜在变量的个体估计没有直接意义，我们应该只关心我们估计整个函数/回归的效果如何。


相反，我们认为检查 BART 变量的收敛性是贝叶斯工作流程的重要组成部分。使用 pmb.plot_convergence 的主要原因是 BART 变量通常是一个大向量（我们估计每个观测值的分布），因此我们需要检查大量诊断。此外，1.01 的 R-hat 阈值不是硬阈值，选择此值是假设检查一个或几个 R-hat（并且链足够长以准确估计其自相关），如果我们观察到大量 R-hat，即使我们的推断没有问题，也预计其中一些会大于 1.01 阈值（或我们选择的任何阈值）。因此，公平的分析应包括多重比较调整，而这正是 pmb.plot_convergence 为您自动执行的操作。那么，如何读取其输出？我们有两个面板，一个用于 ESS，一个用于 R-hat。蓝线是这些值的经验累积分布，对于 ESS，我们希望整条曲线都在虚线上方，对于 R-hat，我们希望曲线完全在虚线下方。在上图中，我们可以看到我们勉强达到了 ESS，并且对于 R-hat，我们只有极少数值高于阈值。我们的结果没用了吗？很可能不是。但为了确定，我们可能需要进行更多抽样。






偏依赖图#


为了帮助我们解释模型的结果，我们将使用偏依赖图。这是一种显示一个协变量对预测变量的边际影响的图。也就是说，协变量 \(X_i\) 对 \(Y\) 的影响是什么，同时我们对所有其他协变量 (\(X_j, \forall j \not = i\)) 求平均值。这种类型的图并非 BART 独有。但它们经常在 BART 文献中使用。PyMC-BART 提供了一个实用函数，可以从 BART 随机变量生成此图。






pmb.plot_pdp(μ, X=X, Y=Y, grid=(2, 2), func=np.exp, var_discrete=[3]);








../_images/f33e0ab5a969a86584bc4ebb18e3b6888779cfa5990d9766abeaf09b1e0e7740.png





从这张图中，我们可以看到每个协变量对预测值的主要影响。这非常有用，因为我们可以恢复超出单调递增或递减效应的复杂关系。例如，对于 hour 协变量，我们可以看到 8 点和 17 点左右有两个峰值，午夜时分有一个最小值。


在解释偏依赖图时，我们应该注意此图中的假设。首先，我们假设变量是独立的。例如，当计算 hour 的影响时，我们必须边缘化 temperature 的影响，这意味着为了计算 hour=0 时的偏依赖值，我们包括所有观察到的温度值，这可能包括午夜时未观察到的温度，因为较低的温度比较高的温度更可能出现。我们只看到平均值，因此如果对于一个协变量，一半的值与预测变量正相关，另一半与预测变量负相关。偏依赖图将是平坦的，因为它们的贡献将相互抵消。这是一个可以通过使用个体条件期望图 pmb.plot_ice(...) 来解决的问题。请注意，所有这些假设都是偏依赖图的假设，而不是我们模型的假设！事实上，BART 可以轻松适应变量的交互作用（尽管 BART 中的先验正则化了高阶交互作用）。有关解释机器学习模型的更多信息，您可以查看“可解释的机器学习”一书 [Molnar, 2019]。


最后，与其他回归方法一样，我们应该注意我们在单个变量上看到的效果是以包含其他变量为条件的。因此，例如，虽然 humidity 看起来大多是平坦的，这意味着此协变量对使用的自行车数量影响很小。这可能是因为 humidity 和 temperature 在某种程度上是相关的，一旦我们在模型中包含 temperature，humidity 就不会提供太多额外信息。例如，尝试再次拟合模型，但这次以 humidity 作为单个协变量，然后再次拟合模型，以 hour 作为单个协变量。您应该看到此单变量模型的结果与先前图中 hour 协变量的结果非常相似，但与 humidity 协变量的结果不太相似。






变量重要性#


正如我们在上一节中看到的那样，偏依赖图可以可视化地让我们了解每个协变量对预测结果的贡献程度。此外，PyMC-BART 提供了一种新颖的方法来评估模型中每个变量的重要性。您可以在下图中看到一个示例。


在 x 轴上，我们有协变量的数量，在 y 轴上，我们有完整模型（包含所有变量）和受限模型（仅包含变量子集）的预测值之间的 R²（皮尔逊相关系数的平方）。


在本例中，最重要的变量是 hour，其次是 temperature、humidity，最后是 workingday。请注意，第一个 R² 值是仅包含变量 hour 的模型的值，第二个 R² 值是包含两个变量（hour 和 temperature）的模型的值，依此类推。除了此排名之外，我们可以看到，即使是具有单个组件 hour 的模型也非常接近完整模型。更重要的是，具有两个组件 hour 和 temperature 的模型在平均意义上与完整模型没有区别。误差条表示来自后验预测分布的 94% HDI。这意味着我们应该期望仅包含 hour 和 temperature 的模型与包含四个变量（hour、temperature、humidity 和 workingday）的模型具有相似的预测性能。






vi_results = pmb.compute_variable_importance(idata_bikes, μ, X)
pmb.plot_variable_importance(vi_results);








../_images/fba22ea65fe421a05029893c942f4e1062ce62ca6713e6ddc8bdfabff59d7ba4.png





plot_variable_importance 速度很快，因为它做出了两个假设


    

  • 变量的排名是通过一个简单的启发式方法计算的。我们只是计算一个变量在所有回归树中被包含的次数。直觉是，如果一个变量很重要，那么它在拟合树中出现的频率应该高于不太重要的变量。
    • 

  • 用于计算 R² 的预测来自已拟合的树。例如，为了估计具有变量 hour 的 BART 模型的效果，我们修剪不包含此变量的分支。这使得计算速度更快，因为我们不需要找到一组新的树。
    • 



除了使用“计数启发式”之外。它还可以执行向后搜索，pmb.plot_variable_importance(..., method="backward")。在内部，这将计算完整模型的 R²，然后计算所有变量少于完整模型的模型的 R²，然后计算所有变量少于两个的模型，依此类推。在每个阶段，我们都会丢弃给出最低 R² 的变量。向后方法会更慢，因为我们需要计算更多模型的预测。






样本外预测#


在本节中，我们想要展示如何使用 BART 进行样本外预测。我们将使用与之前相同的数据集，但这次我们将数据拆分为训练集和测试集。我们将使用训练集来拟合模型，并使用测试集来评估模型。




回归#


让我们首先将此数据建模为回归问题。在这种情况下，我们将数据随机拆分为训练集和测试集。






X_train, X_test, Y_train, Y_test = train_test_split(X, Y, test_size=0.2, random_state=RANDOM_SEED)









现在，我们拟合与上述相同的模型，但这次使用协变量的共享变量，以便我们可以替换它们以生成样本外预测。






with pm.Model() as model_oos_regression:
    X = pm.MutableData("X", X_train)
    Y = Y_train
    α = pm.Exponential("α", 1)
    μ = pmb.BART("μ", X, np.log(Y))
    y = pm.NegativeBinomial("y", mu=pm.math.exp(μ), alpha=α, observed=Y, shape=μ.shape)
    idata_oos_regression = pm.sample(random_seed=RANDOM_SEED)
    posterior_predictive_oos_regression_train = pm.sample_posterior_predictive(
        trace=idata_oos_regression, random_seed=RANDOM_SEED
    )









/home/osvaldo/anaconda3/envs/pymc/lib/python3.11/site-packages/pymc/data.py:321: FutureWarning: MutableData is deprecated. All Data variables are now mutable. Use Data instead.
  warnings.warn(
Multiprocess sampling (4 chains in 4 jobs)
CompoundStep
>NUTS: [α]
>PGBART: [μ]









Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 40 seconds.
The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details
Sampling: [y]













接下来，我们替换模型中的数据并从后验预测分布中采样。






with model_oos_regression:
    X.set_value(X_test)
    posterior_predictive_oos_regression_test = pm.sample_posterior_predictive(
        trace=idata_oos_regression, random_seed=RANDOM_SEED
    )









Sampling: [y, μ]













最后，我们可以将后验预测分布与观察到的数据进行比较。






 隐藏代码单元格来源



fig, ax = plt.subplots(
    nrows=2, ncols=1, figsize=(8, 7), sharex=True, sharey=True, layout="constrained"
)

az.plot_ppc(
    data=posterior_predictive_oos_regression_train, kind="cumulative", observed_rug=True, ax=ax[0]
)
ax[0].set(title="Posterior Predictive Check (train)", xlim=(0, 1_000))

az.plot_ppc(
    data=posterior_predictive_oos_regression_test, kind="cumulative", observed_rug=True, ax=ax[1]
)
ax[1].set(title="Posterior Predictive Check (test)", xlim=(0, 1_000));










../_images/01973d5c4c2d57e5cd65e5394cc5cad470f8cbb43217f648bfad08b3b5c5ef77.png





耶！结果看起来相当合理 🙂！






时间序列#


我们可以使用 hour 特征从时间序列的角度查看相同的数据。从这个角度来看，我们需要确保我们不打乱数据，以免泄漏信息。因此，我们使用 hour 特征定义训练集-测试集拆分。






train_test_hour_split = 19

train_bikes = bikes.query("hour <= @train_test_hour_split")
test_bikes = bikes.query("hour > @train_test_hour_split")

X_train = train_bikes[features]
Y_train = train_bikes["count"]

X_test = test_bikes[features]
Y_test = test_bikes["count"]









然后，我们可以运行相同的模型（但使用不同的输入数据！）并生成如上所述的样本外预测。






with pm.Model() as model_oos_ts:
    X = pm.MutableData("X", X_train)
    Y = Y_train
    α = pm.Exponential("α", 1 / 10)
    μ = pmb.BART("μ", X, np.log(Y))
    y = pm.NegativeBinomial("y", mu=pm.math.exp(μ), alpha=α, observed=Y, shape=μ.shape)
    idata_oos_ts = pm.sample(random_seed=RANDOM_SEED)
    posterior_predictive_oos_ts_train = pm.sample_posterior_predictive(
        trace=idata_oos_ts, random_seed=RANDOM_SEED
    )









/home/osvaldo/anaconda3/envs/pymc/lib/python3.11/site-packages/pymc/data.py:321: FutureWarning: MutableData is deprecated. All Data variables are now mutable. Use Data instead.
  warnings.warn(
Multiprocess sampling (4 chains in 4 jobs)
CompoundStep
>NUTS: [α]
>PGBART: [μ]









Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 52 seconds.
The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details
Sampling: [y]













我们生成样本外预测。






with model_oos_ts:
    X.set_value(X_test)
    posterior_predictive_oos_ts_test = pm.sample_posterior_predictive(
        trace=idata_oos_ts, random_seed=RANDOM_SEED
    )









Sampling: [y, μ]













与上面类似，我们可以将后验预测分布与观察到的数据进行比较。






 隐藏代码单元格来源



fig, ax = plt.subplots(
    nrows=2, ncols=1, figsize=(8, 7), sharex=True, sharey=True, layout="constrained"
)

az.plot_ppc(data=posterior_predictive_oos_ts_train, kind="cumulative", observed_rug=True, ax=ax[0])
ax[0].set(title="Posterior Predictive Check (train)", xlim=(0, 1_000))

az.plot_ppc(data=posterior_predictive_oos_ts_test, kind="cumulative", observed_rug=True, ax=ax[1])
ax[1].set(title="Posterior Predictive Check (test)", xlim=(0, 1_000));











/home/osvaldo/anaconda3/envs/pymc/lib/python3.11/site-packages/IPython/core/events.py:93: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.
  func(*args, **kwargs)
/home/osvaldo/anaconda3/envs/pymc/lib/python3.11/site-packages/IPython/core/pylabtools.py:152: UserWarning: Creating legend with loc="best" can be slow with large amounts of data.
  fig.canvas.print_figure(bytes_io, **kw)




../_images/ccf03f5803ab380b6309c1786838946f16478178807cccb4950ddf0bccdd3e5f.png





哇！这看起来不太对劲！测试集上的预测看起来很奇怪 🤔。为了更好地理解发生了什么，我们可以将预测绘制为时间序列






 隐藏代码单元格来源



fig, ax = plt.subplots(figsize=(12, 6))
az.plot_hdi(
    x=X_train.index,
    y=posterior_predictive_oos_ts_train.posterior_predictive["y"],
    hdi_prob=0.94,
    color="C0",
    fill_kwargs={"alpha": 0.2, "label": r"94$\%$ HDI (train)"},
    smooth=False,
    ax=ax,
)
az.plot_hdi(
    x=X_train.index,
    y=posterior_predictive_oos_ts_train.posterior_predictive["y"],
    hdi_prob=0.5,
    color="C0",
    fill_kwargs={"alpha": 0.4, "label": r"50$\%$ HDI (train)"},
    smooth=False,
    ax=ax,
)
ax.plot(X_train.index, Y_train, label="train (observed)")
az.plot_hdi(
    x=X_test.index,
    y=posterior_predictive_oos_ts_test.posterior_predictive["y"],
    hdi_prob=0.94,
    color="C1",
    fill_kwargs={"alpha": 0.2, "label": r"94$\%$ HDI (test)"},
    smooth=False,
    ax=ax,
)
az.plot_hdi(
    x=X_test.index,
    y=posterior_predictive_oos_ts_test.posterior_predictive["y"],
    hdi_prob=0.5,
    color="C1",
    fill_kwargs={"alpha": 0.4, "label": r"50$\%$ HDI (test)"},
    smooth=False,
    ax=ax,
)
ax.plot(X_test.index, Y_test, label="test (observed)")
ax.axvline(X_train.shape[0], color="k", linestyle="--", label="train/test split")
ax.legend(loc="center left", bbox_to_anchor=(1, 0.5))
ax.set(
    title="BART model predictions for bike rentals",
    xlabel="observation index",
    ylabel="number of rentals",
);










../_images/775b81cc935aac704876e8223e50d8c7dda74de63f43548c20bdedf90fd76272.png





此图帮助我们理解测试集上性能不佳的原因：回想一下，在初始模型的变量重要性排名中，我们看到 hour 是最重要的预测变量。另一方面，我们的训练数据只看到 hour 值直到 \(19\)（因为这是我们的训练集-测试集阈值）。由于 BART 学习如何划分（训练）数据，因此它无法区分例如 \(20\) 和 \(22\) 之间的 hour 值。它只关心这两个值都大于 \(19\)。这在使用 BART 时非常重要！这解释了为什么在存在趋势成分的情况下不应将 BART 用于时间序列预测。在这种情况下，最好先对数据进行去趋势化处理，用 BART 对剩余部分进行建模，并用不同的模型对趋势进行建模。










作者#


    

  • 由 Osvaldo Martin 于 2021 年 12 月编写 (pymc-examples#259)
    • 

  • 由 Osvaldo Martin 于 2022 年 5 月更新 (pymc-examples#323)
    • 

  • 由 Osvaldo Martin 于 2022 年 9 月更新
    • 

  • 由 Osvaldo Martin 于 2022 年 11 月更新
    • 

  • Juan Orduz 于 2023 年 1 月添加了样本外部分
    • 

  • 由 Osvaldo Martin 于 2023 年 3 月更新
    • 

  • 由 Osvaldo Martin 于 2023 年 11 月更新
    • 

  • 由 Osvaldo Martin 于 2024 年 12 月更新
    • 







参考文献#







[1]

Miriana Quiroga, Pablo G Garay, Juan M. Alonso, Juan Martin Loyola, and Osvaldo A Martin. Bayesian additive regression trees for probabilistic programming. 2022. URL: https://arxiv.org/abs/2206.03619, doi:10.48550/ARXIV.2206.03619.





[2]

Christoph Molnar. Interpretable Machine Learning. Christoph Molnar, 2019. URL: https://christophm.github.io/interpretable-ml-book/.





[3]

Osvaldo A Martin, Ravin Kumar, and Junpeng Lao. Bayesian Modeling and Computation in Python. Chapman and Hall/CRC, 2021. doi:10.1201/9781003019169.












水印#






%load_ext watermark
%watermark -n -u -v -iv -w









Last updated: Mon Dec 23 2024

Python implementation: CPython
Python version       : 3.11.5
IPython version      : 8.16.1

pymc_bart : 0.6.0
pymc      : 5.19.1
arviz     : 0.20.0.dev0
pandas    : 2.1.2
matplotlib: 3.8.4
numpy     : 1.26.4

Watermark: 2.4.3













许可声明#


此示例库中的所有笔记本均根据 MIT 许可证 提供，该许可证允许修改和再分发用于任何用途，前提是版权和许可声明得以保留。






引用 PyMC 示例#


要引用此笔记本，请使用 Zenodo 为 pymc-examples 存储库提供的 DOI。




重要提示


许多笔记本都改编自其他来源：博客、书籍…… 在这种情况下，您也应该引用原始来源。


另请记住引用您的代码使用的相关库。




这是一个 bibtex 中的引用模板


@incollection{citekey,
  author    = "<notebook authors, see above>",
  title     = "<notebook title>",
  editor    = "PyMC Team",
  booktitle = "PyMC examples",
  doi       = "10.5281/zenodo.5654871"
}





渲染后可能看起来像