Python 贝叶斯分析第三版（二）

原文：annas-archive.org/md5/23c902bea72dd81eede00e0b57b70813

译者：飞龙

协议：CC BY-NC-SA 4.0

第六章

使用 Bambi 建模

一件好工具能改善你的工作方式，一件伟大的工具能改善你的思维方式。——Jeff Duntemann

在第四章中，我们描述了线性回归模型的基本成分以及如何将其推广以更好地适应我们的需求。在这一章中，我们将继续学习线性模型，但这次，我们将使用 Bambi [Capretto et al., 2022]，一个基于 PyMC 构建的高层贝叶斯模型构建接口。Bambi 旨在使拟合线性模型，包括分层模型，变得极其简单。我们将看到 Bambi 的领域不仅限于线性模型。

我们将学习以下内容：

• 使用 Bambi 构建并拟合模型

• 使用 Bambi 分析结果

• 多项式回归与样条函数

• 分布模型

• 分类预测变量

• 交互作用

• 使用 Kulprit 进行变量选择

6.1 一种语法解决所有问题

PyMC 有一个非常简单而富有表现力的语法，它允许我们构建任意模型。这通常是一个福音，但有时也可能成为负担。Bambi 则专注于回归模型，这种限制使得语法和功能更加专注，正如我们将看到的那样。

Bambi 使用了类似许多 R 包（如 nlme、lme4 和 brms）所用的 Wilkinson 公式语法。假设data是一个像表格 6.1中展示的那样的 pandas DataFrame。

	y	x	z	g
0	-0.633494	-0.196436	-0.355148	A 组
1	2.32684	0.0163941	-1.22847	B 组
2	0.999604	0.107602	-0.391528	C 组
3	-0.119111	0.804268	0.967253	A 组
4	2.07504	0.991417	0.590832	B 组
5	-0.412135	0.691132	-2.13044	C 组

表格 6.1：一个示例 pandas DataFrame

使用这些数据，我们想要构建一个从x预测y的线性模型。使用 PyMC，我们会做类似以下代码块中的模型：

代码 6.1

with pm.Model() as lm: 
    Intercept = pm.Normal("Intercept", 0, 1) 
    x = pm.Normal("x", 0, 1) 
    y_sigma = pm.HalfNormal("sigma", 1) 
    y_mean = Intercept + x * data["x"] 
    y = pm.Normal("y", y_mean, y_sigma, observed=data["y"])

Bambi 使用的公式语法让我们可以更紧凑地定义一个等效的模型：

代码 6.2

a_model = bmb.Model("y ∼ x", data)

在波浪符号（∼）的左侧是因变量，右侧是自变量（们）。通过这种语法，我们只是指定了均值（在 PyMC 模型lm中的μ）。默认情况下，Bambi 假设似然函数是高斯分布；你可以通过family参数来更改这一点。公式语法并没有指定先验分布，只是描述了因变量和自变量之间的关系。Bambi 会自动为我们定义（非常）弱的信息先验。我们可以通过打印 Bambi 模型来获取更多信息。如果你打印a_model，你应该会看到类似如下的内容：

    Formula: y ~ x      Family: gaussian          Link: mu = identity  Observations: 117      Priors:  target = mu      Common-level effects          Intercept ~ Normal(mu: 0.02, sigma: 2.8414)          x ~ Normal(mu: 0.0, sigma: 3.1104)        Auxiliary parameters          sigma ~ HalfStudentT(nu: 4.0, sigma: 1.1348)

第一行显示了我们用来定义模型的公式，第二行是似然函数。第三行是链接函数。然后我们有用于拟合模型的观测数量，接下来告诉我们我们正在线性建模高斯的参数mu。输出的后半部分显示了模型结构：共同水平效应，在本例中是截距(Intercept)和斜率(x)，以及辅助参数，即所有非线性建模的参数，本例中是高斯标准差。

您可以通过将字典传递给bmb.Model的priors参数来覆盖默认先验分布。例如，如果我们想为变量x的系数和辅助参数sigma定义自定义先验分布，我们可以这样做：

代码 6.3

priors = {"x": bmb.Prior("HalfNormal", sigma=3), 
          "sigma": bmb.Prior("Gamma",  mu=1, sigma=2), 
          } 
a_model_wcp = bmb.Model("y ∼ x", data, priors=priors)

结果，我们将得到以下模型规范：

       Formula: y ~ x          Family: gaussian            Link: mu = identity    Observations: 117          Priors:      target = mu          Common-level effects              Intercept ~ Normal(mu: 0.02, sigma: 2.837)              x ~ HalfNormal(sigma: 3.0)            Auxiliary parameters              sigma ~ Gamma(mu: 1.0, sigma: 2.0)

如果您想从模型中省略截距，可以这样做：

代码 6.4：

no_intercept_model = bmb.Model("y ∼ 0 + x", data)

或者甚至这样做：

代码 6.5：

no_intercept_model = bmb.Model("y ∼ -1 + x", data)

打印模型no_intercept_model，您会看到截距不再存在。

如果我们想要包含更多变量怎么办？我们可以这样做：

代码 6.6

model_2 = bmb.Model("y ∼ x + z", data)

我们还可以包含组级效应（层次结构）；例如，如果我们想要使用变量g来部分汇总x的估计值，我们可以这样做：

代码 6.7

model_h = bmb.Model("y ∼ x + z + (x | g)", data)

我们可以在图 6.1中看到这个模型的视觉表示。注意变量1|g_offset和x|g_offset。默认情况下，Bambi 拟合一个非居中的分层模型；您可以通过参数noncentered来更改这一点。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file169.png

图 6.1：model_h的视觉表示

公式语法非常简单，但也非常强大。我们只是初步展示了可以使用它做什么。如果您想深入了解，可以查看公式文档 bambinos.github.io/formulae/。 formulae 是负责解析威尔金森公式用于 Bambi 的 Python 包。

6.2 自行车模型，Bambi 版本

我们将首先使用自行车模型来说明如何使用 Bambi，该模型来自第 4章。我们可以通过以下方式加载数据：

代码 6.8

bikes = pd.read_csv("data/bikes.csv")

现在我们可以构建和拟合模型：

代码 6.9

model_t = bmb.Model("rented ∼ temperature", bikes, family="negativebinomial") 
idata_t = model_t.fit()

图 6.2展示了模型的视觉表示。如果您想要直观地检查先验分布，可以使用model.plot_priors()：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file170.png

图 6.2：用命令model.graph()计算的自行车模型的视觉表示

现在让我们绘制后验均值和后验预测分布（预测）。为了使图表看起来漂亮，省略了一些细节，做这些操作的代码是：

代码 6.10

_, axes = plt.subplots(1, 2, sharey=True, figsize=(12, 4)) 
bmb.interpret.plot_predictions(model_t, idata_t, 
                               "temperature", ax=axes[0]) 
bmb.interpret.plot_predictions(model_t, idata_t, 
                               "temperature", pps=True, ax=axes[1])

plot_predictions 是 Bambi 子模块 interpret 中的一个函数。这个函数通过绘制条件调整预测，帮助分析回归模型，展示（条件）响应分布中的某个参数如何随（某些）插值的解释变量变化。我们可以在图 6.3中看到这段代码的结果。左侧面板显示后验均值和 94% HDI，而右侧面板显示后验预测分布（租用自行车的预测分布）。请注意，预测的不确定性远大于均值的不确定性（pps=False）。这是因为后验预测分布考虑了模型参数的不确定性以及数据的不确定性，而均值的后验分布仅考虑了截距和斜率参数的不确定性。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file171.png

图 6.3：自行车模型的后验均值和后验预测分布

当我们有多个解释变量时，plot_cap 的实用性变得更加明显。例如，假设我们拟合一个模型，使用温度和湿度来预测租赁的自行车数量：

代码 6.11

model_th = bmb.Model("rented ∼ temperature + humidity", bikes, 
                     family="negativebinomial") 

idata_th = model_th.fit() 

bmb.interpret.plot_predictions(model_th, idata_th, ["temperature", "humidity"], 
                               subplot_kwargs={"group":None, "panel":"humidity"})

在图 6.4中，我们可以看到五个面板，每个面板展示了在不同湿度值下，租赁自行车数量随温度变化的情况。正如你所看到的，租赁的自行车数量随温度上升而增加，但在湿度较低时，斜率更大。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file172.png

图 6.4：带有温度和湿度的自行车模型的后验均值

6.3 多项式回归

使用线性回归模型拟合曲线的一种方法是构建一个多项式，像这样：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file173.jpg

我们称m为多项式的次数。

有两点需要注意。首先，多项式回归仍然是线性回归；线性是指系数（β），而不是变量（x）。第二点需要注意的是，我们正在凭空创建新变量。唯一被观测到的变量是x，其余的只是x的幂。通过观察到的变量创造新变量是回归中完全有效的“技巧”；有时这种变换可以通过理论来解释或说明（比如取婴儿身长的平方根），但有时它仅仅是为了拟合曲线。多项式的直观理解是，对于给定的x值，阶数越高的多项式，曲线的灵活性越强。1 阶多项式是直线，2 阶多项式是可以向上或向下的曲线，3 阶多项式是可以先上升再下降的曲线（或反之），以此类推。注意我说的是“可以”，因为如果我们有一个 3 阶多项式，比如 β[0] + β[1]x + β[2]x² + β[3]x³，但系数β[2]和β[3]是 0（或几乎为 0），那么曲线就会是一条直线。

使用 Bambi 定义多项式回归有两种方式。我们可以写出原始多项式：

代码 6.12

"y ∼ x + I(x ** 2) + I(x ** 3) + I(x ** 4)"

在这里，我们使用恒等函数I()来明确表示我们希望将x提升到某个幂次。我们需要这样做，因为**操作符在 Bambi 中有特殊的含义。如果我们使用这种语法，实际上是在告诉 Bambi 将y的均值建模为 α + β[0]x + β[0]x² + β[0]x³ + β[0]x⁴。

或者，我们可以写：

代码 6.13

"y ∼ poly(x, 4)"

这也会生成一个 4 阶的多项式，但多项式项之间是正交的，这意味着项之间的相关性被降低了。不深入数学细节，这至少有两个重要的结果与标准多项式相比。首先，估计可能在数值上更稳定；其次，系数的解释方式不同。在标准多项式回归中，系数可能很难解释，因为改变一个系数的值会影响整个多项式。相比之下，正交多项式可以让你更清晰地解释每个项的影响，因为它们彼此独立。尽管系数的解释不同，但其他结果保持不变。例如，使用这两种方法你应该得到相同的预测结果。

让我们构建一个 4 阶的正交多项式来拟合自行车数据。在这个例子中，我们将使用hour变量：

代码 6.14

model_poly4 = bmb.Model("rented ∼ poly(temperature, degree=4)", bikes, 
                       family="negativebinomial") 
idata_poly4 = model_poly4.fit()

图 6.5 显示了后验均值和后验预测分布。在第一行，你将看到一个 1 阶的多项式，它相当于一个线性模型。在第二行，你将看到一个 4 阶的多项式。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file174.png

图 6.5：温度和湿度下的自行车模型的后验均值和后验预测分布

多项式的一个问题是它们作用全局。当我们应用一个阶数为m的多项式时，我们实际上是在说自变量与因变量之间的关系在整个数据集中都是阶数m。当数据的不同区域需要不同灵活性时，这可能会造成问题。例如，这可能导致曲线过于灵活。随着阶数的增加，拟合变得更加敏感于数据点的移除，或者等同于对未来数据的加入。换句话说，随着阶数的增加，模型更容易发生过拟合。贝叶斯多项式回归通常较少遭遇这种“过度”灵活性，因为我们通常不使用平坦的先验，并且我们不会计算单一的系数集，而是计算整个后验分布。不过，我们还是可以做得更好。

6.4 样条曲线

一般来说，写出非常灵活的模型的方法是将函数B[m]应用于X[m]，然后将它们与系数β[m]相乘：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file175.jpg

我们可以自由选择B[m]，例如，我们可以选择多项式。但我们也可以选择其他函数。一种常见的选择是使用 B 样条；我们不打算讨论它们的定义，但可以将它们视为一种创建平滑曲线的方式，使得我们在获得像多项式一样的灵活性时，过拟合的可能性较小。我们通过使用分段多项式来实现这一点，即将多项式限制在仅影响数据的一部分的区域内。图 6.6 显示了逐步增加阶数的分段多项式的三个示例。虚线垂直线表示“节点”，它们是用于限制区域的点，虚灰线表示我们要逼近的函数，黑线表示分段多项式。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file176.png

图 6.6：逐步增加阶数的分段多项式

图 6.7 显示了 1 阶和 3 阶样条的示例；底部的点表示“节点”，虚线表示 B 样条。在顶部，我们展示了所有具有相等权重的 B 样条；我们使用灰度来突出显示我们有多个 B 样条。在底部面板中，每个 B 样条的权重不同（我们通过β[m]系数对它们进行加权）；如果我们对加权后的 B 样条进行求和，最终得到的黑线就是结果。这条黑线通常被称为“样条曲线”。我们可以使用贝叶斯统计来找到 B 样条的适当权重。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file177.png

图 6.7：1 阶（分段线性）或 3 阶（立方样条）的 B 样条及其结果样条曲线。

我们可以通过使用bs函数，在 Bambi 中使用 B 样条。例如，让我们对自行车数据拟合一个 3 阶的样条曲线：

代码 6.15

num_knots = 6 
knots = np.linspace(0, 23, num_knots+2)[1:-1] 
model_spline = bmb.Model("rented ∼ bs(hour, degree=3, knots=knots)", bikes, 
                         family="negativebinomial") 
idata_spline = model_spline.fit()

图 6.8 显示了租赁自行车数量在深夜时最低。然后有一个增加，可能是因为人们醒来去上班、上学或做其他活动。我们在大约第 8 小时左右有一个第一个高峰，然后略微下降，接着在大约第 18 小时出现第二个高峰，可能是因为人们下班回家，之后出现稳定的下降。注意，曲线并不是非常平滑；这并不是因为样条模型的原因，而是因为数据的原因。我们的测量是在离散的时间点（每小时）进行的。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file178.png

图 6.8：样条模型的后验均值

在处理样条时，我们必须做出一个重要的决策，那就是确定结点的数量和位置。这可能是一个相当令人头疼的任务，因为最佳的结点数量及其间距并不立刻显现出来。一个有用的建议是考虑基于分位数而非均匀分布来确定结点的位置——像这样 knots = np.quantile(bikes.hour, np.linspace(0, 1, num_knots))。通过这样做，我们会在数据量较大的区域放置更多的结点，而在数据量较少的区域放置较少的结点。这样可以得到一个更具适应性的近似，从而有效地捕捉数据点密度较高区域的变化。此外，我们可能希望尝试使用不同数量和位置的结点来拟合样条，然后评估结果，使用诸如 LOO 等工具，就像我们在第五章中看到的那样。

6.5 分布模型

我们之前看到，可以使用线性模型来表示均值（或位置参数）以外的其他参数。例如，我们可以为高斯分布的均值和标准差分别使用线性模型。这些模型通常被称为分布模型。分布模型的语法非常相似；我们只需要为我们想要建模的辅助参数添加一行。例如，高斯分布的σ，或负二项分布的α。

现在让我们重现第四章中的一个例子——婴儿示例：

代码 6.16

formula = bmb.Formula( 
    "length ∼ np.sqrt(month)", 
    "sigma ∼ month" 
) 
model_dis = bmb.Model(formula, babies) 
idata_dis = model_dis.fit()

图 6.9 显示了 model_dis（变化的 sigma）和常数 sigma 模型的后验分布值。我们可以看到，当 sigma 允许变化时，得到的值既低于也高于常数 sigma 的估计值，这意味着当我们不允许 sigma 改变时，我们会低估或高估这个参数。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file179.png

图 6.9：婴儿数据的常数与变化的 sigma

图 6.10 显示了 model_dis 的后验拟合。注意，模型能够捕捉到婴儿成长过程中变化的变异性。这个图与图 4.13非常相似。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file180.png

图 6.10：model_dis的后验拟合

在使用 PyMC 时，我们看到从后验预测分布中进行采样时，需要定义“Xs”作为可变数据，然后在计算后验预测分布之前更新该变量。而在 Bambi 中，这不需要。我们可以使用predict方法，通过将新值传递给data参数来预测新值。例如，让我们预测一只小企鹅在 0.5 个月（15 天）时的体长：

代码 6.17

model_dis.predict(idata_dis, kind="pps", data=pd.DataFrame({"month":[0.5]}))

6.6 类别预测变量

类别变量表示不同的组或类别，这些类别只能从有限的集合中取值。通常，这些值是标签或名称，本身不具备数值意义。以下是一些例子：

• 政治立场：保守派、自由派或进步派。

• 性别：女性或男性。

• 客户满意度：非常不满意、不满意、中立、满意或非常满意。

线性回归模型可以轻松地处理类别变量；我们只需要将类别编码为数字。可以通过几种方式做到这一点，Bambi 可以轻松地为我们处理这些细节。真正的难点在于结果的解释，我们将在接下来的两节中深入探讨。

6.6.1 类别企鹅

对于当前的例子，我们将使用 palmerpenguins 数据集，Horst 等人[2020]，该数据集包含 344 个观测值和 8 个变量。目前，我们关心的是将企鹅的体重建模为喙长的函数。预计随着喙长的增加，企鹅的体重也会增加。本例的新颖之处在于，我们将考虑类别变量species。在这个数据集中，物种变量有 3 个类别或级别，分别是 Adelie、Chinstrap 和 Gentoo。图 6.11 显示了我们要建模的变量的散点图。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file181.png

图 6.11：3 种企鹅的喙长与体重的关系

让我们加载数据并拟合模型：

代码 6.18

penguins = pd.read_csv("data/penguins.csv").dropna() 

model_p = bmb.Model("body_mass ∼ bill_length + species", data=penguins) 
idata_p = model_p.fit()

请注意，定义 Bambi 模型中的类别变量时没有特殊的语法。Bambi 可以自动检测并处理它们。

图 6.12 显示了model_p的森林图。注意到什么意外情况了吗？Adelie 没有后验值。这不是错误。默认情况下，Bambi 将具有 N 个级别（3 个物种）的类别变量编码为 N-1 个虚拟变量（2 个物种）。因此，物种-Chinstrap 和物种-Gentoo 的系数被建模为与基线模型的偏差：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file182.jpg

为了更清楚地说明这一点，我们来查看几个图表。我们可以通过图 6.12 理解为，Chinstrap 的体重平均比 Adelie 轻 0.89。Gentoo 也是如此，不过这次我们需要在基线模型的均值上加上 0.66。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file183.png

图 6.12: 来自 model_p 的森林图

你可以通过查看图 6.13 来验证这两个陈述的真实性。注意中间是阿德利（Adelie），下方是拟企鹅（Chinstrap，-0.89），上方是根趾企鹅（Gentoo，0.58），它们的三条线基本平行。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file184.png

图 6.13: 来自 model_p 的平均样本预测

6.6.2 关于分层模型的关系

在第 3 章中，我们讨论了并对比了汇集模型和分层（或部分汇集）模型。我们在那里展示了通常情况下，我们利用数据的结构或层次。遵循该章节的逻辑，你可以认为阿德利、根趾企鹅和拟企鹅虽然是不同的物种，但都是企鹅。因此，分层建模它们的体重可能是个好主意。你这样想是正确的。那么这种模型与我们在本节中使用的模型有什么区别呢？

区分因素在于斜率和截距的微妙组成部分。在后者的情况下，斜率在三种企鹅物种中保持不变，而截距可以变化：阿德利的Intercept + 0，拟企鹅的Intercept + species[Chinstrap]，根趾企鹅的Intercept + species[Gentoo]。因此，该模型突出了不同的截距，同时保持了统一的斜率。

如果我们建立了层次模型 body_mass ~(bill_length|species)，我们将请求部分汇集的斜率和截距。如果我们建立了模型 body_mass ~(0 + bill_length | species)，我们将请求部分汇集的斜率和公共截距。

除了这些特定的模型外，在考虑将预测变量作为分组变量或分类预测变量使用时，通常有必要问问该变量是否包括所有可能的类别（例如所有星期的所有天，所有物种等），还是仅包括子组（某些学校或少数音乐流派）。如果我们有所有可能的类别，那么我们可能更喜欢将其建模为分类预测变量，否则建模为分组变量。

正如我们之前讨论过的，通常在决定哪一个模型最佳之前，我们会创建多个模型。最佳模型是与你分析目标一致、提供有意义的见解，并准确代表数据潜在模式的模型。探索多个模型，使用适当的标准（例如第 5 章中讨论的标准）比较它们的性能，并考虑每个模型对研究或决策过程的实际影响通常是个好主意。

6.7 交互作用

交互效应或统计交互效应发生在独立变量对响应的影响因另一个独立变量的值而变化时。交互作用可以发生在两个或多个变量之间。一些例子包括:

• 教育水平和收入的影响：较高的教育水平可能对某一性别的收入产生更强的正向影响，从而形成教育和性别之间的交互作用。

• 药物疗效与年龄：一种对年长者效果更好的药物，而对年轻人效果较差。

• 运动和饮食对减重的影响：对于做很少或不做运动的人来说，饮食对减重的影响可能较小；而对于做中等强度运动的人来说，饮食对减重的影响可能较大。

• 作物生长的温度和湿度：一些作物可能在高温高湿的环境中茁壮成长，而其他作物则可能在较凉爽、湿度较低的环境中表现更好。

当两个或更多变量的联合效应不等于它们单独效应的总和时，我们就有了交互作用。因此，如果我们有如下模型，我们就无法建模交互作用：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file185.jpg

建模交互作用效应的最常见方法是将两个（或更多）变量相乘。例如，考虑以下这样的模型：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file186.jpg

在建模交互作用效应时，通常还会包括主效应/项。

新的预测因子

将两个变量相乘可以看作是一种技巧，类似于我们在多项式回归（或对给定变量的任何变换）中使用的方法。我们不是将一个预测因子与自身相乘，而是将两个不同的预测因子相乘，得到一个新的预测因子。

在 PyMC 模型中定义两个变量之间的交互作用非常简单；我们只需要将这两个预测因子相乘，并添加一个系数。对于 Bambi 模型来说，更加简单；我们使用 : 运算符。为了让区别更加清晰，让我们看看一个有交互作用和没有交互作用的模型示例：

代码 6.19

# No interaction 
model_noint = bmb.Model("body_mass ∼ bill_depth + bill_length", 
                        data=penguins) 

#Interaction 
model_int = bmb.Model("body_mass ∼ bill_depth + bill_length + 
                      bill_depth:bill_length", 
                      data=penguins) 

idata_noint = model_noint.fit() 
idata_int = model_int.fit()

我们现在使用 Bambi 的 plot_prediction 来比较不同 bill_length 值如何作为 bill_depth 的函数影响 body_mass。图 6.14 显示了结果。我们为 bill_depth 在 5 个固定的 bill_length 值下计算了平均回归拟合。左侧是 model_noint（无交互作用）结果，右侧是 model_int（有交互作用）结果。我们可以看到，当没有交互作用时，bill_depth 的拟合曲线在不同的 bill_length 水平上是平行的。相反，当存在交互作用时，这些曲线不再平行，正因为 bill_depth 的变化对 body_mass 变化的影响不再是常数，而是受到 bill_length 值的调节。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file187.png

图 6.14：model_noint（左）和 model_int（右）的样本内平均预测

如果你生成了像图 6.14 这样的图，但不是固定 bill_length，而是决定固定 bill_depth，你将观察到类似的行为。

在本书的 GitHub 仓库中（github.com/aloctavodia/BAP3），你将找到文件 interactions.ipynb。这个脚本生成一个 3D 图形，我希望它能帮助你建立关于添加交互项时我们所做工作的直觉。如果你运行它，你会看到，当没有交互项时，我们拟合的是一个 2D 平面，一张像纸一样的平坦表面。但是，当添加交互项时，你会拟合一个弯曲的表面。将 interactions.ipynb 的结果与 图 6.14 进行比较。

我们刚刚通过可视化看到，解释包含交互项的线性模型并不像解释没有交互项的线性模型那样简单。让我们从数学上来看一下。

假设我们有一个包含 2 个变量 X[0] 和 X[1] 以及它们之间的交互作用的模型：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file188.jpg

我们可以将这个模型重写为：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file189.jpg

或者甚至像这样：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file190.jpg

从这个表达式中，我们可以看到：

• 交互项可以理解为线性模型中的线性模型。

• 交互作用是对称的；我们可以将其理解为 X[0] 作为 X[1] 的函数的斜率，同时也是 X[1] 作为 X[0] 的函数的斜率。这也可以从交互式图形中看到。

• 我们之前知道，β[0] 系数可以解释为 μ 对 X[0] 单位变化的响应量（这就是我们称之为斜率的原因）。如果我们添加一个交互项，那么只有当 X[1] = 0 时，这个关系才成立。试着使用交互式图形自己看看。数学上，这是真的，因为当 X[1] = 0 时，β[2]X[1] = 0，因此 X[0] 的斜率简化为 β[0]X[0]。通过对称性，同样的推理可以应用于 β[1]。

6.8 使用 Bambi 解释模型

在本章中，我们已经多次使用了 bmb.interpret_plot_predictions。但这并不是 Bambi 提供的唯一帮助我们理解模型的工具。另一个工具是 bmb.interpret_plot_comparisons。这个工具帮助我们回答这样的问题：“当我们比较某个变量的两个值时，保持其他所有值不变，预测差异是什么？”

我们使用上一节的 model_int，因此无需重新拟合模型。我们使用以下代码块生成 图 6.15：

代码 6.20

bmb.interpret.plot_comparisons(model_int, idata_int, 
                               contrast={"bill_depth":[1.4, 1.8]}, 
                               conditional={"bill_length":[3.5, 4.5, 5.5]})

图 6.15 显示了当将假设的企鹅的 bill_depth 从 1.8 比较到 1.4 时，期望的差异是：

• 大约 0.8 千克，当喙长为 3.5 厘米时

• -0.6 千克，当喙长为 4.5 厘米时

• 大约 -2 千克，当喙长为 5.5 厘米时

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file191.png

图 6.15：在 3 个固定的 bill_length 值下，将 bill_depth 从 1.8 厘米对比到 1.4 厘米

如果你希望以表格形式查看信息，可以使用函数bmb.interpret.comparisons，这样你将得到一个 DataFrame 而不是图表。

另一个有用的函数是bmb.interpret_plot_slopes，它可用于计算给定值下的“瞬时变化率”或斜率。我们使用以下代码块生成图 6.16：

代码 6.21

bmb.interpret.plot_slopes(model_int, idata_int, 
                          wrt={"bill_depth":1.8}, 
                          conditional={"bill_length":[3.5, 4.5, 5.5]},

图 6.16 显示了bill_depth为 1.8 时的斜率：

• ≈ 2 kg/cm，适用于鸟嘴长度为 3.5 cm 时

• -1.4 kg/cm，适用于鸟嘴长度为 4.5 cm 时

• ≈ -5 kg/cm，适用于鸟嘴长度为 5.5 cm 时

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file192.png

图 6.16：对于bill_length的三个固定值，bill_depth在 1.8 cm 时的斜率

如果你希望以表格形式查看信息，可以使用函数bmb.interpret.slopes，这样你将得到一个 DataFrame 而不是图表。

在本节中，我们仅仅触及了bmb.interpret模块工具的表面。这个模块是 Bambi 的一个非常有用的功能，尤其适用于含有交互作用和/或具有除恒等函数之外链接函数的模型。我强烈建议你阅读 Bambi 文档，获取更多未在此介绍的示例和细节。

6.9 变量选择

变量选择是指从一组潜在预测变量中识别出对模型结果最相关的变量的过程。我们进行变量选择时假设，只有一部分变量对感兴趣的结果有显著影响，而其他变量几乎没有或没有额外的价值。

构建模型时，最“贝叶斯”的做法可能是将我们可能想到的所有变量都包含在一个模型中，然后利用该模型的后验分布进行预测或了解变量之间的关系。这是最“贝叶斯”的方法，因为我们尽可能地使用数据，并在后验中纳入变量重要性的未知性。然而，比贝叶斯更“贝叶斯”并不总是最佳选择。我们在第五章中已经看到，即使贝叶斯因子是贝叶斯定理的直接结果，它们也可能存在问题。

当以下情况出现时，执行变量选择是一个好主意：

• 我们需要降低测量成本。例如，在医学中，我们可能有资金和资源进行一项试点研究，测量 200 名患者的 30 个变量。但是我们无法对数千人做同样的事情。或者，我们可能能在开阔地带放置大量传感器来更好地建模作物收益，但我们无法将其扩展到一个国家的规模。降低成本并不总是意味着金钱或时间；在与人类或其他动物工作时，减少疼痛和不适也很重要。例如，我们可能想预测患者发生心脏病发作的风险。我们可以通过测量很多变量来做到这一点，但我们也可以通过测量少量更不具侵入性的变量来实现。

• 我们希望减少计算成本。这对于小型和简单模型来说不是问题，但当我们有大量变量、大量数据或两者兼有时，计算成本可能会变得无法承受。

• 我们寻求更好地理解重要的关联结构。也就是说，我们有兴趣理解哪些变量提供更好的预测。需要明确的是，我们不是在讨论因果关系。虽然统计模型，尤其是广义线性模型（GLMS），可以用来推断因果关系，但这需要额外的步骤和假设。在本书中，我们不讨论如何进行因果推断。如果你想要一个非常基础的因果推断介绍，请参见这个视频：www.youtube.com/watch?v=gV6wzTk3o1U。如果你更加严肃，您可以查阅 Scott Cunningham 的在线书籍《因果推断：混音带》[Cunningham，2021] mixtape.scunning.com/。

• 当我们希望获得一个对数据生成分布变化更具韧性的模型时，我们可以将变量选择视为一种使模型对不具代表性数据更加稳健的方法。

6.9.1 投影预测推断

有许多方法可以执行变量选择。在本节中，我们将重点介绍其中一种叫做投影预测推断的方法[Piironen 等人，2020，McLatchie 等人，2023]。我们专注于这种方法的主要原因是，它在广泛的领域中表现出了非常好的性能。

投影预测推断的主要步骤如下：

1. 生成一个参考模型，即包含所有你认为可能相关和/或你能够测量的变量的模型。

2. 生成一组子模型，即仅包含参考模型中某些子集变量的模型。

3. 将参考模型的后验分布投影到子模型中。

4. 选择一个预测结果与参考模型足够接近的最小模型。

在进行投影预测推断时，我们只需要执行一次贝叶斯推断，仅针对参考模型。对于子模型，后验分布是通过投影得到的。不深入讨论技术细节，投影过程是通过某种方式找到子模型的参数，使得子模型的预测尽可能接近参考模型的预测。这个投影过程可以以计算效率高的方式进行，因此估计后验分布的成本比使用 MCMC 方法低几个数量级。这一点很重要，因为当我们增加参考模型中的变量数时，可能的子模型总数会呈指数增长。考虑一下我们需要评估所有可能的组合，不重复变量。例如，假设我们有四个变量（A，B，C 和 D），需要评估 7 个模型，分别是 A，B，C，AB，BC，AC 和参考模型 ABC。七个模型听起来不多，但当我们增加到 8 个变量时，我们将需要评估 92 个不同的模型。看到没有，我们将变量数翻倍，模型的数量却增加了 10 倍以上！

当然，减少需要探索的子模型总数是有办法的。例如，我们可以使用某些廉价方法筛选出最有前景的变量，然后只对这些变量进行投影预测推断。另一个替代方案被称为前向搜索；也就是说，我们首先拟合与我们所拥有的变量数量相等的模型。然后选择一个模型/变量，即生成与参考模型的预测最接近的那个模型。接着，我们生成所有包含上一阶段选择变量的两个变量子模型，依此类推。如果我们对一个包含 8 个变量的参考模型进行这种前向程序，而不是 92 个不同的模型，我们只需要评估 36 个模型。

另一个在进行投影预测推断时需要考虑的因素是，我们只为参考模型提供了先验分布。子模型没有明确的先验分布；它们只是通过投影过程继承了参考模型的先验分布。

投影预测在实践中有效的原因之一，正是得益于参考模型的使用。通过将子模型拟合到参考模型在样本中的预测，而不是拟合观察到的数据，我们能够过滤掉数据中的噪声。这有助于将更相关的变量与不太相关的变量区分开。另一个因素是，在选择子模型时使用了交叉验证，如在第五章中所讨论的那样。

6.9.2 使用 Kulprit 进行投影预测

Kulprit 是一个用于投影预测推断的 Python 包。它与 Bambi 配合使用，我们可以传递一个用 Bambi 构建的参考模型，Kulprit 会为我们完成所有复杂的工作。为了说明如何使用 Kulprit，我们将使用体脂数据集 [Penrose et al., 1985]。这个数据集包含了 251 个个体的测量数据，包括他们的年龄、体重、身高、腹围等。我们的目的是预测体脂百分比（通过 siri 变量估算）。由于获得准确的体脂测量既昂贵又可能对患者造成困扰，我们希望在减少测量数量的同时保持 siri 的良好预测准确性。原始数据集包含 13 个变量；为了简化示例，我已经预先选择了 6 个变量。

我们首先需要像往常一样定义并拟合一个 Bambi 模型。我们必须确保包含参数 idata_kwargs=’log_likelihood’:True。Kulprit 内部会计算 ELPD，正如我们在 第五章 中讨论的那样，我们需要在 InferenceData 对象中包含对数似然值，以便能够估算 ELPD：

代码 6.22

model = bmb.Model("siri ∼ age + weight + height + abdomen + thigh + wrist", 
                  data=body) 
idata = model.fit(idata_kwargs={'log_likelihood': True})

完成这些后，我们就可以开始使用 Kulprit 了。首先，我们需要调用 ProjectionPredictive 类，并传入 Bambi 模型和从该模型拟合得到的 idata。然后我们让 Kulprit 执行搜索；默认情况下，它会进行前向搜索：

代码 6.23：

ppi = kpt.ProjectionPredictive(model, idata) 
ppi.search()

搜索完成后，我们可以让 Kulprit 根据 ELPD 比较各个子模型。子模型会按 ELPD 从低到高排序，如 图 6.17 所示。在 x 轴上，我们有子模型的大小，即变量的数量；我们从零开始，因为我们包含了仅有截距的模型。虚线灰色线表示参考模型的 ELPD。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file193.png

图 6.17：使用 Kulprit 得到的子模型比较。由 ppi.plot_compare 生成。

然后我们可以看到，大小为 3 的子模型几乎与参考模型等效。但具体包含了哪些变量呢？如果我们在执行搜索后打印 ppi 对象，我们将得到一个子模型公式的有序列表，列表顺序与通过命令 ppi.plot_compare 得到的图中的顺序相匹配：

代码 6.24

print(ppi)

  0 siri ~ 1
  1 siri ~ abdomen
  2 siri ~ abdomen + wrist
  3 siri ~ abdomen + wrist + height
  4 siri ~ abdomen + wrist + height + age
  5 siri ~ abdomen + wrist + height + age + weight
  6 siri ~ abdomen + wrist + height + age + weight + thigh

然后我们可以看到，大小为 3 的模型是包含变量 abdomen（腹部）、wrist（手腕）和 height（身高）的模型。这个结果告诉我们，如果我们想选择一个比参考模型更少变量的模型，但预测精度相似，那么这是一个不错的选择。根据不同的情境，其他子模型也可能是一个好选择。例如，我们可能会认为大小为 2 和 3 的子模型之间的差异非常小。因此，我们可能愿意牺牲一些精度来选择一个更小的模型。对于这个例子来说，测量患者身高可能不会带来太大问题，但在其他场景下，添加第三个变量可能会很昂贵、麻烦、危险等等。

另一种解释 图 6.17 的方式是注意到大小为 3 或更大的模型的 ELPD 值是非常接近的。可能的情况是，如果我们使用稍微不同的数据集，甚至是相同的数据集，但增加更多的后验样本，我们可能会得到一个略有不同的顺序。因此，如果我们有许多大小为 3 的模型，它们可能具有相同的实际预测精度，我们可以通过外部因素来为选择第三个变量提供理由，比如它的测量是否容易或便宜，或者哪个对患者来说更不痛苦等等。总之，和其他统计工具一样，结果不应盲目接受，而是要结合上下文来解读；你应该有最终的决定权，工具应帮助你做出决策。

好的，假设我们确实对 Kulprit 计算的大小为 3 的子模型感兴趣；我们可以通过以下方式获取它：

代码 6.25

submodel = ppi.project(3)

从 submodel 对象中，我们可以检索一些有用的信息，比如 Bambi 的模型 submodel.model 或 InferenceData 对象 submodel.idata。

对于解释这两个对象，有一点需要注意——submodel.model 是一个由公式生成的 Bambi 模型。因此，它的先验将是 Bambi 自动计算的那些。但 Kulprit 计算的后验，存储在 submodel.idata.posterior 中，并非直接来自该模型。相反，它是使用投影预测推理（而非 MCMC）计算的，先验是在投影步骤中隐式继承的（而非显式先验）。图 6.18 展示了这种投影后的后验。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file194.png

图 6.18：大小为 3 的子模型的投影后验

我们可以信任预测的后验分布吗？在非常一般的条件下，这应该是一个有效的后验分布，因此我们可以信任它。它应该足够提供参数值的大致概念，当然，这对于变量选择也是足够的。缺乏显式的先验可能会使模型的解释更加困难，但如果你只关心预测，这应该不是问题。当然，你始终可以使用 Bambi（或 PyMC）像往常一样显式地计算完整的后验，并在需要时自己指定先验。图 6.19 显示了通过 Bambi（真实）计算的子模型后验的森林图和通过 Kulprit（预测）近似的子模型后验。请注意，这里有两个可能的差异来源：MCMC 方法与投影预测方法之间的内在差异，以及两个模型的不同先验。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file195.png

图 6.19：通过 Kulprit 计算的子模型（siri ~abdomen + wrist + height）后验与通过 Bambi 计算的参考模型后验的比较；未在两个模型中共享的变量已被省略

Kulprit 是一个非常新的库，将不断发展，用户可以期待不久后会有大量的增强和改进。如果 Kulprit 引起了你的兴趣，你可以通过报告问题、提出建议、改进文档或参与其代码库的开发来帮助其发展，网址是github.com/bambinos/kulprit。

6.10 小结

在本章中，我们已经看到如何使用 Bambi 拟合贝叶斯模型，作为纯 PyMC 模型的替代方案。我们从最简单的情况开始，一个包含单一预测变量的模型，然后过渡到更复杂的模型，包括多项式、样条、分布模型、包含分类预测变量的模型以及交互作用模型。

Bambi 的主要优势在于它非常易于使用；它与 R 的formula语法非常相似。内部，Bambi 定义了弱信息量的先验，并处理了复杂模型中可能繁琐的细节。主要的缺点是它不像 PyMC 那样灵活。Bambi 能够处理的模型范围是 PyMC 模型范围的一个小子集。不过，这个子集包含了工业界和学术界中最常用的许多统计模型。Bambi 的优势不仅仅在于轻松构建模型，还在于更容易的模型解释。在本章中，我们已经看到如何使用 Bambi 的interpret模块更好地理解我们拟合的模型。最后，我们还看到了如何使用 Kulprit 进行投影预测推断并进行变量选择。投影预测推断为变量选择提供了一种有前景的方法，而 Kulprit 则是一种有前景的 Python 化方式。

6.11 练习

1. 阅读 Bambi 文档（bambinos.github.io/bambi/）并学习如何指定自定义先验。

2. 应用前一点学到的内容，为 model_t 的斜率指定一个 HalfNormal 先验分布。

3. 定义一个类似 model_poly4 的模型，但使用 raw 多项式，比较两个模型的系数和均值拟合效果。

4. 用你自己的话解释什么是分布模型。

5. 将 model_spline 扩展为一个分布模型。使用另一个样条来建模 NegativeBinomial 家族的 α 参数。

6. 创建一个名为 model_p2 的模型，用于预测 body_mass，其预测变量为 bill_length、bill_depth、flipper_length 和 species。

7. 使用 LOO 方法比较前一点中的模型和 model_p。

8. 使用 interpret 模块中的函数来解释 model_p2，并同时使用图表和表格。

加入我们的社区 Discord 空间

加入我们的 Discord 社区，与志同道合的人一起学习，并与超过 5000 名成员共同成长，访问链接：packt.link/bayesian

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file1.png

第七章

混合模型

…父亲有狮子的形态，母亲有蚂蚁的形态；父亲吃肉，母亲吃草。于是他们孕育了蚂蚁狮子…——《虚构生物志》

拉普拉塔河（也称为拉普拉塔河或拉普拉塔江）是地球上最宽的河流，也是阿根廷和乌拉圭之间的自然边界。在 19 世纪末，这条河沿岸的港口区是土著人、非洲人（大多数是奴隶）和欧洲移民的混居地。这种文化交汇的一个结果是，欧洲音乐如华尔兹和马祖卡与非洲坎东贝舞曲和阿根廷米隆加（后者又是非洲裔美国节奏的混合）相融合，创造了我们现在所称的探戈舞和音乐。

将先前存在的元素混合在一起是创造新事物的好方法，这不仅仅是在音乐领域。在统计学中，混合模型是建模的常见方法之一。这些模型通过混合更简单的分布来获得更复杂的分布。例如，我们可以结合两个高斯分布来描述双峰分布，或者结合多个高斯分布来描述任意分布。尽管使用高斯分布非常普遍，但原则上我们可以混合任何我们想要的分布族。混合模型用于不同的目的，例如直接建模子人群，或者作为处理无法用简单分布描述的复杂分布的有用技巧。

本章将涵盖以下主题：

• 有限混合模型

• 零膨胀模型和障碍模型

• 无限混合模型

• 连续混合模型

7.1 理解混合模型

当总体人群是不同子人群的组合时，自然会出现混合模型。一个常见的例子是给定成人群体中的身高分布，它可以被描述为女性和男性子人群的混合。另一个经典的例子是手写数字的聚类。在这种情况下，期望有 10 个子人群是非常合理的，至少在十进制系统中是这样！如果我们知道每个观察值属于哪个子人群，通常来说，利用这些信息将每个子人群建模为一个独立的群体是一个好主意。然而，当我们无法直接访问这些信息时，混合模型就显得特别有用。

分布的混合

许多数据集无法用单一的概率分布准确描述，但可以通过将其描述为多种分布的混合来进行建模。假设数据来自混合分布的模型被称为混合模型。

在构建混合模型时，并不一定需要相信我们正在描述数据中的真实子群体。混合模型也可以作为一种统计技巧，为我们的工具箱增添灵活性。以高斯分布为例，我们可以将其作为许多单峰且近似对称分布的合理近似。那么，对于多峰或偏斜分布呢？我们能用高斯分布来建模吗？当然可以，如果我们使用高斯混合。

在高斯混合模型中，每个成分都是一个均值不同的高斯分布，并且通常（但不一定）具有不同的标准差。通过组合高斯分布，我们可以为我们的模型增添灵活性，并拟合复杂的数据分布。事实上，我们可以通过适当组合高斯分布来逼近几乎任何我们想要的分布。分布的具体数量将取决于近似的准确度和数据的细节。实际上，我们在本书中的许多图表中已经使用了这个思想。核密度估计（KDE）技术就是这一思想的非贝叶斯实现。从概念上讲，当我们调用 az.plot_kde 时，函数会在每个数据点上方放置一个固定方差的高斯分布，然后将所有单独的高斯分布求和，以逼近数据的经验分布。图 7.1 显示了如何将 8 个高斯分布混合以表示复杂分布的示例，就像一条蟒蛇消化一只大象，或者一个帽子，取决于你的视角。

在 图 7.1 中，所有高斯分布具有相同的方差，并且它们都集中在灰色的点上，这些点代表可能的未知人群中的样本点。如果你仔细观察，可能会发现两个高斯分布重叠在一起。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file196.png

图 7.1：作为高斯混合模型的 KDE 示例

无论我们是相信子群体的存在，还是将其作为数学上的便利（甚至是介于两者之间的某种观点），混合模型通过使用分布的混合来描述数据，是为我们的模型增添灵活性的一种有效方式。

7.2 有限混合模型

构建混合模型的一种方法是考虑两个或更多分布的有限加权混合。然后，观察数据的概率密度是 K 个子群体的概率密度加权求和：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file197.jpg

我们可以将 w[i] 解释为成分 i 的概率，因此其值被限制在区间 [0, 1] 内，并且它们的和需要为 1。成分 p(y|θ[i]) 通常是简单的分布，比如高斯分布或泊松分布。如果 K 是有限的，那么我们就得到了一个有限混合模型。为了拟合这样的模型，我们需要提供一个 K 值，无论是因为我们事先知道正确的值，还是因为我们能做出有根据的猜测。

从概念上讲，解决混合模型的问题，我们只需要正确地将每个数据点分配到一个成分。在一个概率模型中，我们可以通过引入一个随机变量来实现这一点，随机变量的功能是指定某个观察值被分配到哪个成分。这个变量通常被称为潜变量，因为我们不能直接观察到它。

对于只有两个成分的混合模型（K = 2），我们可以使用抛硬币问题模型作为构建模块。在该模型中，我们有两种可能的结果，并且使用伯努利分布来描述它们。由于我们不知道正面或反面的概率，我们使用 Beta 分布作为先验分布。如果我们把硬币的正面和反面换成任何两个组（或成分或类别），我们可以在得到 0 时将观察值分配到一个组，得到 1 时分配到另一个组。这一切都很好，但在混合模型中，我们可能有两个或更多组，因此我们需要使这个想法更加通用。我们可以通过注意到伯努利分布到 K 个结果的推广是类别分布，而 Beta 分布到更高维度的推广是狄利克雷分布来实现这一点。如果我们分配的成分是高斯分布，就像在图 7.1 中一样，那么图 7.2 显示了一个类似 Kruschke 风格的模型图。圆角框表示我们有 K 个成分，类别变量决定了我们用哪个成分来描述给定的数据点。请注意，只有 μ[k] 依赖于不同的成分，而 σ[μ] 和 σ[σ] 是所有成分共享的。这只是一个建模选择；如果需要，我们可以改变它，并允许其他参数依赖于每个成分。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file198.png

图 7.2：带有高斯成分的混合模型的 Kruschke 风格图

我们将在 PyMC 中实现这个模型，但在此之前，让我介绍一下类别分布和狄利克雷分布。如果你已经熟悉这些分布，可以跳过接下来的两节，直接跳到示例部分。

7.2.1 类别分布

类别分布是最一般的离散分布，它由一个向量进行参数化，其中每个元素指定每个可能结果的概率。图 7.3 表示了类别分布的两个可能实例。点表示类别分布的值，而连续的虚线是一个视觉辅助，帮助我们轻松理解分布的形状：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file199.png

图 7.3：类别分布族的两个成员。虚线只是视觉辅助。

7.2.2 狄利克雷分布

Dirichlet 分布存在于单纯形中，你可以将其视为一个 n 维三角形；一个 1-单纯形是一个线段，2-单纯形是一个三角形，3-单纯形是一个四面体，依此类推。为什么是单纯形？直观上，因为这个分布的输出是一个K维向量，其元素被限制在区间 [0,1] 内，并且它们的和为 1。正如我们所说，Dirichlet 分布是 Beta 分布的推广。因此，理解 Dirichlet 分布的一个好方法是将其与 Beta 分布进行比较。Beta 分布用于处理两种结果的问题：一种的概率是p，另一种的概率是 1 − p。正如我们所见，p + (1 − p) = 1。Beta 分布返回一个两元素向量（p, q = 1 − p），但在实际应用中，我们省略了q，因为一旦知道p的值，结果就完全确定了。如果我们要将 Beta 分布扩展到三种结果，我们需要一个三元素向量（p, q, r），其中每个元素都在区间 [0,1] 内，且它们的和为 1。类似于 Beta 分布，我们可以使用三个标量来对这种分布进行参数化，我们可以称它们为α, β 和 γ；然而，由于希腊字母只有 24 个，我们很容易就会用完它们。相反，我们可以使用一个名为α的K维向量。请注意，我们可以将 Beta 和 Dirichlet 看作是比例分布。图 7.4展示了当k = 3 时，Dirichlet 分布的 4 个成员。在顶部是概率密度函数（pdf），在底部是来自该分布的样本。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file200.png

图 7.4：Dirichlet 分布族的四个成员

好的，现在我们已经具备了实现我们第一个混合模型的所有组件。

7.2.3 化学混合物

为了让事情更具体，我们来做一个例子。我们将使用我们在第 2 章中看到的化学位移数据。图 7.5展示了这组数据的直方图。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file201.png

图 7.5：化学位移数据的直方图

我们可以看到，这些数据无法用像高斯这样的单一分布来准确描述，但如果我们使用多个分布，就能做到这一点，就像我们在图 7.1中展示的那样；也许三个或四个分布就能解决问题。虽然有充分的理论依据（我们在这里不讨论），表明这些化学位移数据来自于 40 个子群体的混合分布，但仅通过观察数据，我们似乎无法恢复出真正的群体，因为它们之间有大量的重叠。

以下代码块展示了在 PyMC 中实现两个分量的高斯混合模型：

代码 7.1

K = 2 
with pm.Model() as model_kg: 
   p = pm.Dirichlet('p', a=np.ones(K)) 
   z = pm.Categorical('z', p=p, shape=len(cs_exp)) 
   means = pm.Normal('means', mu=cs_exp.mean(), sigma=10, shape=K) 
   sd = pm.HalfNormal('sd', sigma=10) 

   y = pm.Normal('y', mu=means[z], sigma=sd, observed=cs_exp) 
   idata_kg = pm.sample()

如果你运行这段代码，你会发现它非常慢，且轨迹看起来非常糟糕（参考第十章了解更多诊断信息）。我们能让这个模型跑得更快吗？可以，让我们看看如何操作。

在model_kg中，我们已经明确将潜在变量z包含在模型中。对这个离散变量进行采样通常会导致后验采样不良。解决这个问题的一种方法是对模型进行重新参数化，使得z不再是模型的显式组成部分。这种重新参数化方法叫做边缘化（marginalization）。将离散变量边缘化通常能够提高速度并改进采样。不幸的是，这需要一些数学技能，而不是每个人都具备。不过幸运的是，我们不需要自己动手，因为 PyMC 已经包含了一个NormalMixture分布。所以，我们可以将混合模型写成如下形式：

代码 7.2

with pm.Model() as model_mg: 
    p = pm.Dirichlet('p', a=np.ones(K)) 
    means = pm.Normal('means', mu=cs_exp.mean(), sigma=10, shape=K) 
    sd = pm.HalfNormal('sd', sigma=5) 
    y = pm.NormalMixture('y', w=p, mu=means, sigma=sd, observed=cs_exp) 
    idata_mg = pm.sample()

让我们通过森林图来检查结果。图 7.6显示了一些有趣的现象。在进入下一个章节之前，花些时间思考一下。你能发现问题吗？我们将在下一节讨论它。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file202.png

图 7.6：model_mg的均值森林图

7.3 混合模型的非可识别性

means参数的形状为 2，从图 7.6我们可以看到，其中一个值大约是 47，另一个接近 57.5。搞笑的是，我们有一条链说means[0]是 47，其他三条链说它是 57.5，means[1]则正好相反。因此，如果我们计算means[0]的均值，我们会得到一个接近 55 的值（57.5 × 3 + 47 × 1），这显然不是正确的值。我们看到的正是一个名为参数非可识别性（parameter non-identifiability）现象的例子。这是因为，从模型的角度来看，如果组件 1 的均值为 47 而组件 2 的均值为 57.5，或者反之，两种情况是等价的。在混合模型的背景下，这也被称为标签交换问题（label-switching problem）。

非可识别性

如果模型的一个或多个参数无法被唯一确定，则该统计模型是不可识别的。如果对于模型参数的多个选择，得到相同的似然函数，则该模型的参数没有被识别。可能出现的情况是数据中没有足够的信息来估计这些参数。在其他情况下，参数可能无法识别，因为模型在结构上不可识别，这意味着即使所有必要的数据都可用，参数仍然无法唯一确定。

对于混合模型，有至少两种方式可以对模型进行参数化，从而消除非可识别性问题。我们可以强制使组件按顺序排列；例如，将组件的均值按严格递增的顺序排列和/或使用信息量大的先验分布。

使用 PyMC，我们可以通过下一个代码块中的转换实现第一个选项。请注意，我们还提供了均值的初始值；任何能确保第一个均值小于第二个均值的方式都可以。

代码 7.3

with pm.Model() as model_mgo: 
    p = pm.Dirichlet('p', a=np.ones(K)) 
    means = pm.Normal('means', mu=cs_exp.mean(), sigma=10, shape=K, 
                      transform=pm.distributions.transforms.ordered, 
                      initval=np.array([cs_exp.mean()-1, cs_exp.mean()+1])) 
    sd = pm.HalfNormal('sd', sigma=10) 
    y = pm.NormalMixture('y', w=p, mu=means, sigma=sd, observed=cs_exp) 

    idata_mgo = pm.sample()

让我们通过森林图检查新的结果。图 7.7 确认我们已解决了不可识别性问题：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file203.png

图 7.7：model_mgo 的均值森林图

7.4 如何选择 K

有关有限混合模型的一个主要问题是如何确定成分的数量。一个经验法则是从相对较小的成分数开始，然后增加成分数来改善模型拟合度。如我们在 第五章 中已经知道的，模型拟合可以通过后验预测检查、ELPD 等指标以及模型者的专业知识来评估。

让我们比较 K = {2*,3,4,5} 的模型。为此，我们将拟合该模型四次，然后保存数据和模型对象以备后用：

代码 7.4

Ks = [2, 3, 4, 5] 

models = [] 
idatas = [] 
for k in Ks: 
    with pm.Model() as model: 
        p = pm.Dirichlet('p', a=np.ones(k)) 
        means = pm.Normal('means', 
                          mu=np.linspace(cs_exp.min(), cs_exp.max(), k), 
                          sigma=cs_exp.var() / k, shape=k, 
                          transform=pm.distributions.transforms.ordered, 
                          ) 
        sd = pm.HalfNormal('sd', sigma=5) 
        y = pm.NormalMixture('y', w=p, mu=means, sigma=sd, observed=cs_exp) 
        idata = pm.sample(random_seed=123, 
                          idata_kwargs={"log_likelihood":True} 
                         ) 

        idatas.append(idata) 
        models.append(model)

图 7.8 展示了 K 个高斯分布的混合模型。黑色实线表示后验均值，灰色线表示后验样本。均值高斯分量使用黑色虚线表示。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file204.png

图 7.8：不同数量高斯分布的高斯混合模型 (K)

从视觉效果来看，K = 2 太低了，但我们如何选择一个更好的值呢？正如我们在 第五章 中讨论过的，我们可以使用后验预测检查来评估感兴趣的测试量，并计算贝叶斯 p 值。图 7.9 展示了此类计算和可视化的示例。K = 5 是最佳解，K = 4 接近最佳解。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file205.png

图 7.9：后验预测检查以选择 K

为了补充后验预测检查，我们可以计算使用 LOO 方法近似的 ELPD。这在 图 7.10 中展示。我们比较的是相同的模型，但 K 值不同。我们可以看到 K = 5 是最佳解，而 K = 4 也接近最佳解。这与 图 7.9 中显示的贝叶斯 p 值一致。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file206.png

图 7.10：使用 LOO 进行模型选择以选择 K

化学位移的例子虽然简单，但展示了有限混合模型的主要思想。在这个例子中，我们使用了高斯分布，因为它们提供了一个良好的数据拟合近似。然而，如果需要，我们也可以使用非高斯分量。例如，我们可以使用：

• 泊松混合模型：假设你在监控每小时进入商店的顾客数量。泊松混合模型可以帮助识别顾客流量的不同模式，例如高峰时段或高峰日，通过假设数据遵循泊松分布的混合。

• 指数混合模型：假设你在研究某种类型的灯泡寿命。指数混合模型可以帮助识别不同寿命的灯泡群体，提示制造质量或环境因素的潜在差异。

在接下来的部分，我们将探索一种非常特殊的混合模型，这种模型涉及两个过程：一个生成零，另一个生成零或非零。

7.5 零膨胀和障碍模型

在计数事物时，比如道路上的汽车、天空中的星星、皮肤上的痣，或几乎任何其他事物，一个选项是不计数某个事物，也就是说得到零。零这个数字通常会由于很多原因出现；比如我们在数红色的车，但是一辆红色的车没有经过街道，或者我们错过了它。如果我们使用泊松分布或负二项分布来建模这种数据，我们会注意到模型生成的零比实际数据少。我们该如何解决这个问题呢？我们可以尝试解决模型预测零比实际观察到的少的具体原因，并将这一因素纳入模型中。但通常情况下，假设我们有两个过程的混合模型可能就足够了，而且更简单：

• 一个过程由一个离散分布建模，概率为 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/Phi_01.png

• 另一个过程以概率 1 − https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/Phi_01.png 生成额外的零

在某些文本中，你会发现 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/Phi_01.png 代表额外的零，而不是 1 − https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/Phi_01.png。这不是什么大问题；只需注意哪个是哪个，尤其是具体示例时。

允许生成“额外”零的分布族被称为零膨胀分布。该家族中最常见的成员包括：

• 零膨胀泊松分布

• 零膨胀负二项分布

• 零膨胀二项分布

在下一节中，我们将使用零膨胀泊松分布来解决回归问题。一旦你了解如何使用这个分布，使用零膨胀负二项分布或零膨胀二项分布将变得非常简单。

7.5.1 零膨胀泊松回归

为了举例说明零膨胀泊松回归模型，我们将使用来自数字研究与教育研究所的数据集（www.ats.ucla.edu/stat/data）。我们有 250 组游客数据，数据包括：每组捕捉到的鱼的数量（count）、每组中的儿童数量（child）以及他们是否带了露营车到公园（camper）。使用这些数据，我们将建立一个模型，通过儿童和露营车变量来预测捕捉到的鱼的数量。

使用 PyMC，我们可以为这个数据编写一个模型，如下所示：

代码 7.5

with pm.Model() as ZIP_reg: 
    <https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/Phi_02.png> = pm.Beta('<https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/Phi_02.png>', 1, 1) 
    *α* = pm.Normal('*α*', 0, 1) 
    *β* = pm.Normal('*β*', 0, 1, shape=2) 
    *θ* = pm.math.exp(*α* + *β*[0] * fish_data['child'] + *β*[1] * fish_data['camper']) 
    yl = pm.ZeroInflatedPoisson('yl', <https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/Phi_02.png>, *θ*, observed=fish_data['count']) 
    trace_ZIP_reg = pm.sample()

camper 是一个二元变量，值为 0 表示不是露营车，1 表示是露营车。表示某个属性是否存在的变量通常被称为虚拟变量或指示变量。请注意，当 camper 的值为 0 时，涉及 β[1] 的项也会变为 0，模型就会简化为一个只有一个自变量的回归模型。我们在第六章 6 中讨论了分类预测变量时已经提到过这一点。

结果如图 7.11所示。我们可以看到，孩子数量越多，捕获的鱼的数量越少。而且，带着露营车旅行的人通常会捕到更多的鱼。如果你检查 child 和 camper 的系数，你会发现我们可以得出以下结论：

• 每增加一个孩子，捕获的鱼的期望数量会减少约 ≈ 0*.*4

• 使用露营车露营会使捕获的鱼的期望数量增加约 2

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file207.png

图 7.11：捕获的鱼数作为孩子数量和是否使用露营车的函数

零膨胀模型与障碍模型密切相关，因此在零膨胀模型的概念仍然清晰的情况下学习障碍模型是非常有益的。

7.5.2 障碍模型

在障碍模型中，伯努利概率决定一个计数变量是否为零或大于零。如果大于零，我们认为已跨越 障碍，这些正值的分布是通过截断在零的分布来确定的。

在混合模型的框架下，我们可以将零膨胀模型视为零和其他某些值（可能是零或非零）的混合。而障碍模型则是零和非零的混合。因此，零膨胀模型只能增加 P(x = 0) 的概率，而对于障碍模型，概率可以变得更小或更大。

在 PyMC 和 Bambi 中可用于障碍模型的分布包括：

• 障碍 Poisson 分布

• 障碍负二项分布

• 障碍伽玛分布

• 障碍对数正态分布

为了说明障碍模型，我们将使用螯虾数据集 [Brockmann, 1996]。螯虾成对地来到海滩进行产卵仪式。此外，孤独的雄性也会来到海岸，聚集在已成对的雌雄螯虾周围，争夺受精卵的机会。这些被称为卫星雄性个体，通常会在特定的巢对附近聚集，忽略其他巢对。我们想要建立模型来预测雄性 卫星 的数量。我们怀疑这个数量与雌性螯虾的特征有关。作为预测变量，我们将使用螯甲的 宽度 和 颜色。螯甲是螯虾的坚硬上壳。颜色则使用从 1 到 4 的整数编码，从浅色到深色。

我们将使用 Bambi 来编码并拟合四个模型。这四个模型的主要区别在于我们将使用四种不同的似然函数或分布族，分别是 Poisson 分布、障碍 Poisson 分布、负二项分布和障碍负二项分布。模型展示在下一个代码块中：

代码 7.6

model_crab_p = bmb.Model("satellite ∼ width + C(color)", 
                         family="poisson", data=crab) 
model_crab_hp = bmb.Model("satellite ∼ width + C(color)", 
                          family="hurdle_poisson", data=crab) 
model_crab_nb = bmb.Model("satellite ∼ width + C(color)", 
                          family="negativebinomial", data=crab) 
model_crab_hnb = bmb.Model("satellite ∼ width + C(color)", 
                           family="hurdle_negativebinomial", data=crab)

请注意，我们已经将颜色编码为 C(color)，以向 Bambi 指示应将其视为分类变量，而不是数值型变量。

图 7.12 显示了我们为马蹄铁数据拟合的四个模型的后验预测检查。灰色条形图代表观察数据的频率。点是根据模型计算的预期值。虚线只是一个视觉辅助工具。我们可以看到，NegativeBinomial 比 Poisson 模型有了改善，而 Hurdle 模型比非膨胀模型有所提升。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file208.png

图 7.12：四个模型在马蹄铁数据上的后验预测检查

图 7.13 显示了通过 LOO 计算的 ELPD 模型比较。Hurdle NegativeBinomial 是最佳模型，而 Poisson 模型是最差的。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file209.png

图 7.13：使用 LOO 进行的四个模型在马蹄铁数据上的模型比较

图 7.12 很好，但还有一种替代的表示方式，叫做悬挂根图 [Kleiber and Zeileis, 2016]，它在诊断和处理计数数据模型中的问题（如过度离散和/或零值过多）时特别有用。请参见图 7.14，这是针对马蹄铁数据和我们的四个模型的一个示例。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file210.png

图 7.14：四个模型在马蹄铁数据上的后验预测检查，附带根图

在悬挂的根图中，我们绘制了观察值和预测值的平方根。这是一种快速的方式，用于大致调整不同计数之间的尺度差异。换句话说，它使得即使是低频率的观察值和预期频率也能更容易进行比较。其次，观察数据的条形图是从预期值“悬挂”下来的，而不是像在图 7.12 中那样从零开始“生长”。由于条形图是悬挂的，如果条形图没有达到零（虚线灰线），则表示模型对该计数的预测过高；如果条形图低于零，则表示模型对该计数的预测过低。

让我们总结一下图 7.12 中每个子图的内容：

• Poisson：零值被低估，1 到 4 的计数被高估。6 及以后的大部分计数也被低估。这个模式表明数据中存在过度离散，而 0 的巨大差异表明零值过多。

• NegativeBinomial：我们可以看到，与 Poisson 模型相比，过度离散得到了更好的处理。我们仍然看到零值被低估，计数 1 和 2 被高估，这可能表明零值过多。

• Hurdle Poisson：正如预期的那样，对于一个 Hurdle 模型，我们对零值的拟合非常完美。对于正值，我们仍然会看到一些偏差。

• Hurdle NegativeBinomial：我们可以看到，该模型能够很好地拟合数据，对于大多数计数，偏差非常小。

7.6 混合模型和聚类

聚类或聚类分析是将对象分组的任务，目的是使得同一组中的对象彼此之间比与其他组的对象更接近。这些组被称为簇，接近程度可以通过多种不同的方式计算，例如使用度量，如欧几里得距离。如果我们选择概率方法，那么混合模型作为解决聚类任务的自然候选者出现。

使用概率模型进行聚类通常被称为基于模型的聚类。使用概率模型可以计算每个数据点属于每个簇的概率。这被称为软聚类，而不是硬聚类，在硬聚类中，每个数据点要么属于一个簇，概率为 0 或 1。我们可以通过引入一些规则或边界将软聚类转化为硬聚类。事实上，你可能还记得，这正是我们将逻辑回归转化为分类方法时所做的操作，其中我们使用 0.5 作为默认边界值。对于聚类，合理的选择是将数据点分配给概率最高的簇。

总结来说，当人们谈论聚类时，通常是在谈论将对象分组，而当人们谈论混合模型时，他们是在谈论使用简单分布的混合来建模更复杂的分布，目的是识别子群体或仅仅为了拥有一个更灵活的模型来描述数据。

7.7 非有限混合模型

对于一些问题，例如尝试对手写数字进行聚类，容易确定我们期望在数据中找到的组数。对于其他问题，我们可以有很好的猜测；例如，我们可能知道我们的鸢尾花样本来自一个只有三种鸢尾花生长的区域，因此使用三个成分是一个合理的起点。当我们不确定成分的数量时，我们可以使用模型选择来帮助我们选择组数。然而，对于其他问题，事先选择组数可能是一个缺点，或者我们可能更感兴趣的是直接从数据中估计这个数量。对于这种类型的问题，贝叶斯解决方案与狄利克雷过程相关。

7.7.1 狄利克雷过程

到目前为止我们所看到的所有模型都是参数模型，意味着它们具有固定数量的参数，我们感兴趣的是估计这些参数，如固定数量的聚类。我们也可以有非参数模型。这些模型的更好名称可能是非固定参数模型或具有可变数量参数的模型，但已经有人为我们决定了名称。我们可以将非参数模型看作是具有理论上无限数量参数的模型。在实践中，我们在某种程度上让数据将理论上无限的参数数量减少到有限的数量。由于数据决定了实际的参数数量，非参数模型非常灵活，并且有可能对欠拟合和过拟合具有鲁棒性。在本书中，我们将看到三种此类模型的例子：高斯过程（GP）、贝叶斯加法回归树（BART）和 Dirichlet 过程（DP）。虽然接下来的章节将分别集中讨论 GP 和 BART，但我们当前的重点将是探索 DPs。

由于 Dirichlet 分布是 Beta 分布的 n 维推广，Dirichlet 过程是 Dirichlet 分布的无限维推广。我知道这一开始可能让人感到困惑，所以在继续之前，花点时间重新阅读上一句。

Dirichlet 分布是在概率空间上的一种概率分布，而 DP 是在分布空间上的概率分布。这意味着从 DP 中单次抽样实际上是一个分布。对于有限混合模型，我们使用 Dirichlet 分布来为固定数量的聚类或组分配先验分布。DP 是一种为非固定数量的聚类分配先验分布的方法。我们可以将 DP 看作是从分布的先验分布中进行抽样的一种方法。

在我们进入实际的非参数混合模型之前，让我们花点时间讨论一些 DP 的细节。DP 的正式定义相对晦涩，除非你非常了解概率论，否则不容易理解，因此让我描述一些 DP 的属性，这些属性对于理解其在非有限混合模型中的作用是非常重要的：

• DP 是一种其实现为概率分布的分布。例如，从高斯分布中，你会抽样数值，而从 DP 中，你会抽样分布。

• DP 由基础分布 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/H.PNG 和 α（一个正实数，称为浓度参数）来指定。α 类似于 Dirichlet 分布中的浓度参数。

• https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/H.PNG 是 DP 的期望值。这意味着 DP 会围绕基础分布生成分布。这在某种程度上等同于高斯分布的均值。

• 随着 α 的增加，实现在分布上的集中度越来越低。

• 在实践中，DP 总是生成离散分布。

• 在极限情况下，α → ∞，DP 的实现将等于基础分布，因此如果基础分布是连续的，DP 将生成一个连续分布。因为这个原因，数学家们说从 DP 生成的分布几乎总是离散的。实际上，alpha 是一个有限数值，因此我们总是使用离散分布。

分布的先验

我们可以将 DP 看作是随机分布 f 的先验，其中基础分布 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/H_gray.PNG 是我们预期 f 应该是什么，而浓度参数 α 表示我们对先验猜测的信心程度。

为了使这些特性更具体，我们再来看一下 图 7.3 中的分类分布。我们可以通过指明 x 轴上的位置和 y 轴上的高度来完全指定这个分布。对于分类分布，x 轴上的位置被限制为整数，并且高度之和必须等于 1。我们保持最后一个限制，但放宽前一个限制。为了生成 x 轴上的位置，我们将从基础分布 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/H.PNG 中采样。原则上，它可以是我们想要的任何分布；因此，如果我们选择高斯分布，位置将是实数线上的任意值。相反，如果我们选择 Beta 分布，位置将限制在区间 [0, 1] 内，如果我们选择泊松分布作为基础分布，位置将限制为非负整数 0, 1, 2, ….

到目前为止都很好，但我们如何选择 y 轴上的值呢？我们遵循一个被称为“想象实验”（Gedanken experiment）的过程，称为“断棒过程”。假设我们有一根长度为 1 的棒子，然后我们把它分成两部分（不一定相等）。我们把一部分放一边，接着把另一部分再分成两部分，然后我们就这样不断地做下去，永远做下去。实际上，由于我们不能真正无限地重复这个过程，我们在某个预定义的值 K 处截断它，但这个基本思路至少在实践中是成立的。为了控制断棒过程，我们使用一个参数 α。当我们增大 α 的值时，我们会把棒子分成越来越小的部分。因此，当 α = 0 时，我们不分棒子，而当 α = ∞ 时，我们会把它分成无限多的部分。图 7.15 展示了从 DP 中抽取的四个样本，分别对应四个不同的 α 值。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file211.png

图 7.15：以高斯分布为基础的断棒过程

从图 7.15中我们可以看到，DP 是一个离散分布。当α增大时，我们从初始单位长度的棒中获得较小的块；注意 y 轴刻度的变化。基础分布（在此图中为 Normal(0, 1)）控制位置。随着α的增加，棒逐渐更多地呈现出基础分布的特征。在本章的配套笔记本中，你将找到生成图 7.15的代码。我强烈建议你玩一下这段代码，以便更好地理解 DP 的直觉。

图 7.1 显示了，如果你在每个数据点上放置一个高斯分布，然后将所有高斯分布相加，你可以近似数据的分布。我们可以使用 DP 做类似的事情，但不是在每个数据点上放置一个高斯分布，而是在每个原始单位长度的棒的的位置上放置一个高斯分布。然后我们通过每块棒的长度来加权每个高斯分布。这个过程为非有限高斯混合模型提供了一般的方法。

或者，我们可以用任何其他分布替代高斯分布，这样我们就得到了一个非有限混合模型的一般方法。图 7.16 展示了这种模型的一个示例。我使用了拉普拉斯分布的混合模型，这只是为了强调你绝不局限于仅使用高斯混合模型：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file212.png

图 7.16：使用 DP 的拉普拉斯混合模型

现在我们已经准备好尝试在 PyMC 中实现 DP。让我们首先定义一个与 PyMC 兼容的stick_breaking函数：

代码 7.7

K = 10 

def stick_breaking(*α*, K): 
    *β* = pm.Beta('*β*', 1., *α*, shape=K) 
    w = *β* * pt.concatenate([[1.], pt.extra_ops.cumprod(1\. - *β*)[:-1]]) + 1E-6 
    return w/w.sum()

我们不是固定α（浓度参数）的值，而是为其定义一个先验。一个常见的选择是 Gamma 分布，如下代码块所示：

代码 7.8

with pm.Model() as model_DP: 
    *α* = pm.Gamma('*α*', 2, 1) 
    w = pm.Deterministic('w', stick_breaking(*α*, K)) 
    means = pm.Normal('means', 
                      mu=np.linspace(cs_exp.min(), cs_exp.max(), K), 
                      sigma=5, shape=K, 
                      transform=pm.distributions.transforms.ordered, 
                     ) 

    sd = pm.HalfNormal('sd', sigma=10, shape=K) 
    obs = pm.NormalMixture('obs', w, means, sigma=sd, observed=cs_exp.values) 
    idata = pm.sample()

因为我们是通过截断的 stick-breaking 过程来近似无限 DP，因此重要的是检查截断值（在此例中为K = 10）是否引入了任何偏差。一种简单的方法是计算每个组件的平均权重，对其进行排序，然后绘制其累积和。为了保险起见，我们应该确保至少有一些组件的权重可以忽略不计；否则，我们必须增加截断值。此类图的示例见图 7.17。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file213.png

图 7.17：DP 组件平均权重的有序累积分布

我们可以看到，只有前 7 个组件在某种程度上是重要的。前 7 个组件占总权重的 99.9%以上（在图 7.17中的灰色虚线），因此我们可以确信选择的上限值（K = 10）对于这组数据来说足够大。

图 7.18 显示了使用 DP 模型估计的平均密度（黑线），以及从后验样本中获得的样本（灰线），以反映估计的不确定性。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file214.png

图 7.18：用于化学位移数据的 DP 混合模型

7.8 连续混合

本章的重点是离散混合模型，但我们也可以有连续混合模型。事实上，我们已经知道其中一些。例如，层次模型也可以解释为连续混合模型，其中每个组的参数来自上层的连续分布。为了更具体一点，想象一下对多个组进行线性回归。我们可以假设每个组都有自己的斜率，或者所有组共享相同的斜率。或者，与其将问题框架设置为两个极端的离散选项，层次模型让我们能够有效地建模这两种选项的连续混合。

7.8.1 一些常见分布是混合分布

BetaBinomial 是一种离散分布，通常用于描述在每次试验成功的概率p未知并假设服从参数α和β的 Beta 分布的情况下，进行n次伯努利试验时成功的次数y：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file215.jpg

也就是说，为了找到观察到结果y的概率，我们需要对所有可能的（且连续的）p值进行平均。因此，BetaBinomial 可以被视为一个连续混合模型。如果 BetaBinomial 模型对你来说听起来很熟悉，那是因为你在书的前两章时已经注意到了！这是我们用于掷硬币问题的模型，尽管我们明确使用了 Beta 和 Binomial 分布，而不是使用已经混合的 Beta-Binomial 分布。

类似地，我们有负二项分布，它可以被理解为伽马-泊松混合模型。也就是说，这是泊松分布的混合，其中速率参数服从伽马分布。负二项分布通常用于解决处理计数数据时常遇到的一个问题。这个问题被称为过度离散。假设你使用泊松分布来建模计数数据，然后你意识到数据中的方差超过了模型的方差；使用泊松分布的一个问题是，均值和方差由同一个参数描述。解决过度离散的一种方法是将数据建模为泊松分布的（连续）混合。通过考虑分布的混合，我们的模型具有更大的灵活性，能够更好地适应数据的均值和方差。

另一种分布混合的例子是 Student’s t 分布。我们将该分布介绍为高斯分布的稳健替代。在这种情况下，t 分布是由均值μ且方差未知的高斯分布的混合所产生的，该方差服从 InverseGamma 分布。

7.9 总结

许多问题可以描述为由不同子群体组成的整体群体。当我们知道每个观察值属于哪个子群体时，就可以将每个子群体特定地建模为一个独立的群体。然而，许多时候我们无法直接获得这些信息，因此使用混合模型来建模这些数据可能是合适的。我们可以利用混合模型来尝试捕捉数据中的真实子群体，或者作为一种通用的统计技巧，通过将简单分布组合来建模复杂的分布。

在本章中，我们将混合模型分为三类：有限混合模型、非有限混合模型和连续混合模型。有限混合模型是由两个或多个分布的有限加权混合组成，每个分布或组分代表数据的一个子群体。原则上，组分可以是我们认为有用的任何东西，从简单的分布（如高斯分布或泊松分布）到更复杂的对象（如层次模型或神经网络）。从概念上讲，为了解决混合模型，我们所需要做的就是将每个数据点正确地分配给一个组分。我们可以通过引入潜在变量z来实现这一点。我们为z使用类别分布，这是一种最通用的离散分布，具有狄利克雷先验，它是 Beta 分布的 n 维推广。对离散变量z进行采样可能会出现问题，因此将其边缘化可能更为方便。PyMC 包括一个正态混合分布和一个执行此边缘化的混合分布，这使得使用 PyMC 构建混合模型变得更加简单。

在本章中，我们研究混合模型时遇到的一个常见问题是，这种模型可能会导致标签交换问题，这是一种非可识别性问题。解决非可识别性的一种方法是强制将各个组分进行排序。有限混合模型的一个挑战是如何确定组分的数量。一种解决方案是对一组模型进行模型比较，基于估计的组分数量来选择最合适的模型。这一估计应当尽可能地借助我们对当前问题的理解。另一个选择是尝试从数据中自动估计组分的数量。为此，我们引入了狄利克雷过程的概念，作为狄利克雷分布的无限维版本，借此我们可以构建一个非参数的混合模型。

最后，为了结束本章，我们简要讨论了许多模型（例如 Beta 二项分布（用于硬币抛掷问题）、负二项分布、学生 t 分布，甚至层次模型）如何被解释为连续混合模型。

7.10 习题

1. 从 3 个高斯分布的混合中生成合成数据。请查阅本章附带的 Jupyter 笔记本，了解如何执行此操作。拟合一个具有 2、3 或 4 个组分的有限高斯混合模型。

2. 使用 LOO 比较练习 1 的结果。

3. 阅读并运行以下关于混合模型的 PyMC 文档示例：

   • 边际高斯混合模型：www.pymc.io/projects/examples/en/latest/mixture_models/marginalized_gaussian_mixture_model.html

   • 相关密度回归：www.pymc.io/projects/examples/en/latest/mixture_models/dependent_density_regression.html

4. 使用负二项分布和障碍负二项分布模型重新拟合 fish_data。使用根图比较这两个模型与本章展示的零膨胀泊松模型。

5. 使用狄利克雷过程重复练习 1。

6. 假设你暂时不知道鸢尾花数据集的正确物种/标签，使用混合模型将三个鸢尾花物种进行聚类，选择一个特征（例如花萼的长度）。

7. 重复练习 6，但这次使用两个特征。

加入我们的社区 Discord 空间

加入我们的 Discord 社区，结识志同道合的人，与超过 5000 名成员一起学习，网址：packt.link/bayesian

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file1.png

第八章

高斯过程

孤单吗？你有你自己。你无限的自己。- Rick Sanchez（至少是 C-137 维度中的那个）

在上一章中，我们学习了狄利克雷过程，这是一种狄利克雷分布的无限维推广，可以用来对未知的连续分布设定先验。在本章中，我们将学习高斯过程，这是一种高斯分布的无限维推广，可以用来对未知的函数设定先验。狄利克雷过程和高斯过程都用于贝叶斯统计中，构建灵活的模型，在这些模型中，参数的数量可以随着数据量的增加而增加。

我们将讨论以下主题：

• 作为概率对象的函数

• 核函数

• 具有高斯似然的高斯过程

• 非高斯似然下的高斯过程

• 希尔伯特空间中的高斯过程

8.1 线性模型与非线性数据

在第四章和第六章中，我们学习了如何构建以下形式的模型：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file216.jpg

这里，θ是某个概率分布的参数，例如，高斯分布的均值，二项分布的p参数，泊松分布的速率，等等。我们称https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/phi.png为逆链接函数，https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/phi.png是我们用来潜在地变换数据的其他函数，比如平方根、多项式函数，或者其他形式。

拟合或学习一个贝叶斯模型可以看作是寻找权重β的后验分布，因此这被称为近似函数的权重视角。正如我们在多项式回归和样条回归中看到的那样，通过让https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/phi.png成为一个非线性函数，我们可以将输入映射到特征空间中。我们还看到，通过使用适当阶数的多项式，我们可以完美地拟合任何函数。但除非我们应用某种形式的正则化，例如，使用先验分布，否则这将导致记忆数据的模型，换句话说，模型的泛化能力非常差。我们还提到，样条回归可以像多项式一样灵活，但具有更好的统计属性。现在，我们将讨论高斯过程，它为通过有效地让数据决定函数复杂度来建模任意函数提供了一个有原则的解决方案，同时避免或至少最小化过拟合的可能性。

以下章节从一个非常实际的角度讨论高斯过程；我们避免了涉及几乎所有相关的数学内容。对于更正式的解释，你可以阅读 Rasmussen 和 Williams 的《Gaussian Processes for Machine Learning》[2005]。

8.2 函数建模

我们将通过首先描述一种将函数表示为概率对象的方法来开始讨论高斯过程。我们可以将一个函数 f 视为从输入集合 X 到输出集合 Y 的映射。因此，我们可以写成：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file217.jpg

表示函数的一种非常粗略的方式是列出每个 x[i] 值对应的 y[i] 值，如 表 8.1 所示。你可能还记得这种表示函数的方式，它来自小学阶段。

x	y
0.00	0.46
0.33	2.60
0.67	5.90
1.00	7.91

表 8.1：一种函数的表格表示（某种程度上）

作为一般情况，X 和 Y 的值将位于实数线上；因此，我们可以将函数视为一个（可能）无限且有序的值对列表 (x[i],y[i])。顺序很重要，因为如果我们打乱这些值的顺序，得到的将是不同的函数。

根据这个描述，我们可以数值地表示我们想要的任何特定函数。但是，如果我们想要概率性地表示函数怎么办呢？那么，我们就需要编码一个概率映射。让我解释一下；我们可以让每个值都是一个具有某种关联分布的随机变量。由于使用高斯分布通常很方便，我们可以假设它们服从具有给定均值和方差的高斯分布。这样，我们就不再描述单一的特定函数，而是描述一个分布族。

为了使讨论更加具体，我们将使用一些 Python 代码来构建并绘制两个这样的函数示例：

代码 8.1

np.random.seed(42) 
x = np.linspace(0, 1, 10) 
y = np.random.normal(0, 1, len(x)) 
plt.plot(x, y, 'o-', label='the first one') 
y = np.zeros_like(x) 

for i in range(len(x)): 
    y[i] = np.random.normal(y[i-1], 1) 
plt.plot(x, y, 'o-', label='the second one') 
plt.legend()

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file218.png

图 8.1：从高斯分布中采样的两个虚拟函数

图 8.1 显示了使用来自高斯分布的样本来编码函数并非那么疯狂或愚蠢，所以我们可能走在正确的道路上。然而，用于生成 图 8.1 的方法是有限的，并且不够灵活。

虽然我们期望实际的函数具有某种结构或模式，但我们表示 第一个 函数的方式并没有让我们编码数据点之间的任何关系。实际上，每个点都是完全独立的，因为我们只是从一个共同的一维高斯分布中随机抽取了 10 个独立的样本。对于 第二个 函数，我们引入了一些依赖关系。点 y[i+1] 的均值是 y[i]，因此我们这里有一定的结构。然而，我们将在接下来的内容中看到，捕捉依赖关系有一种更为通用的方法，而不仅仅是连续点之间的依赖。

在继续之前，让我停下来思考一下，为什么我们使用高斯分布而不是其他概率分布。首先，通过将自己限制为只使用高斯分布，我们在指定不同形状的函数时不会失去任何灵活性，因为每个点可能有自己的均值和方差。其次，从数学的角度来看，使用高斯分布是非常方便的。

8.3 多元高斯和函数

在图 8.1中，我们将一个函数表示为从一维高斯分布中采样得到的集合。另一种选择是使用 n 维多元高斯分布，得到一个长度为n的样本向量。实际上，你可以尝试重现图 8.1，但将np.random.normal(0, 1, len(x))替换为np.random.multivariate_normal，均值为np.zeros_like(x)，标准差为np.eye(len(x))。使用多元正态分布的优势在于，我们可以利用协方差矩阵来编码有关函数的信息。例如，通过将协方差矩阵设置为np.eye(len(x))，我们表明我们在评估函数的 10 个点中，每个点的方差为 1。我们还表明它们之间的方差，即协方差为 0。换句话说，它们是独立的。如果我们将这些零替换为其他数值，我们可以得到描述不同故事的协方差。

我希望你已经开始相信使用多元高斯来表示函数是可能的。如果是这样，那么我们只需要找到一个合适的协方差矩阵。这也是下一节的主题。

8.3.1 协方差函数和核函数

在实践中，协方差矩阵是通过称为核函数的函数来指定的。不幸的是，术语“核”是一个多义词，甚至在统计学文献中也是如此。定义核函数的一种简单方法是，任何返回有效协方差矩阵的函数都可以称为核。但这个定义是自我重复的，且不太直观。一个更具概念性和实用性的定义是，核函数定义了输入空间中数据点之间的相似度度量，而这种相似度决定了一个数据点对预测另一个数据点值的影响程度。

有许多有用的核函数，一个流行的核函数是指数二次核：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file219.jpg

在这里，∥X − X′∥²是平方欧几里得距离：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file220.jpg

对于这个核函数，我们可以看到它是一个对称函数，接受两个输入，如果输入相同则返回 0，否则返回正值。因此，我们可以将指数二次核的输出解释为两个输入之间的相似度度量。

乍一看可能不太明显，但指数二次核的公式与高斯分布非常相似。正因为如此，这个核函数也被称为高斯核函数。术语ℓ被称为长度尺度（或带宽或方差），它控制着核函数的宽度。换句话说，它控制着X值在什么尺度下被认为是相似的。

为了更好地理解核函数的作用，我建议你动手尝试它们。例如，定义一个 Python 函数来计算指数二次核函数：

代码 8.2

def exp_quad_kernel(x, knots, ℓ=1): 
    """exponentiated quadratic kernel""" 
    return np.array([np.exp(-(x-k)**2 / (2*ℓ**2)) for k in knots])

图 8.2 展示了在不同输入下 4 × 4 协方差矩阵的样子。我选择的输入非常简单，包含了值 [−1*,0,1,*2]。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file221.png

图 8.2：输入值（左），协方差矩阵（右）

在 图 8.2 的左面板中，我们展示了输入值。这些是 x 轴上的值，我们将点从 0 标记到 3。因此，点 0 取值 -1，点 1 取值 0，以此类推。在右面板中，我们展示了使用指数二次核函数计算得到的协方差矩阵的热力图。颜色越浅，协方差值越大。如你所见，热力图是对称的，且对角线上的值最大。当我们意识到协方差矩阵中每个元素的值与点之间的距离成反比时，这一点是有意义的，对角线是通过将每个数据点与自身进行比较得到的。最小的值出现在点 0 和点 3，因为它们是最远的两个点。

一旦你理解了这个例子，你应该尝试使用其他输入。请参考本章末的练习 1 以及附带的笔记本（github.com/aloctavodia/BAP3）进行进一步练习。

现在我们更好地理解了如何使用核函数生成协方差矩阵，接下来让我们更进一步，使用协方差矩阵来采样函数。如你在 图 8.3 中所见，高斯核函数意味着一系列不同的函数，而参数 ℓ 控制函数的平滑度。ℓ 的值越大，函数越平滑。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file222.png

图 8.3：四个 ℓ 值的高斯核函数的实现（每个 ℓ 值对应两个实现）

你告诉我你的朋友，我就能告诉你你的未来

核函数将数据点沿 x 轴的距离转化为期望函数值（y 轴）的协方差值。因此，数据点在 x 轴上越接近，我们期望它们在 y 轴上的值越相似。

8.4 高斯过程

现在我们已经准备好理解什么是高斯过程（GPs）以及它们在实践中的应用。根据维基百科，高斯过程的一个正式定义如下：

“由时间或空间索引的随机变量的集合，使得这些随机变量的每个有限集合都具有多元正态分布，即它们的每个有限线性组合服从正态分布。”

这个定义可能并不是很有用，至少在你当前的学习阶段不太有用。理解高斯过程的诀窍在于意识到，高斯过程的概念是一个心理（和数学）框架，因为在实际应用中，我们并不需要直接处理这个无限维的数学对象。相反，我们只在有数据的点上评估高斯过程。通过这样做，我们将无限维的高斯过程压缩成一个有限维的多元高斯分布，维度数与数据点的数量相等。从数学上讲，这种压缩是通过对无限未观察到的维度进行边际化来实现的。理论告诉我们，除了我们观察到的点，忽略（实际上是边际化）所有其他点是可以的。它还保证我们总是会得到一个多元高斯分布。因此，我们可以严格地解释图 8.3 为高斯过程的实际样本！

到目前为止，我们专注于多元正态分布的协方差矩阵，并没有讨论均值。在使用高斯过程时，将多元高斯的均值设置为 0 是常见做法，因为高斯过程足够灵活，可以将均值建模得非常好。但请注意，这样做并没有限制。实际上，对于某些问题，你可能希望参数化地建模均值，而将高斯过程用来建模残差。

高斯过程是对函数的先验

高斯过程是对函数的先验分布，使得在你评估函数的每个点时，它会在该点放置一个具有给定均值和方差的高斯分布。在实践中，高斯过程通常是通过使用核函数来构建的，核函数将 x 轴上的距离转化为 y 轴上的相似度。

8.5 高斯过程回归

假设我们可以将一个变量Y建模为X的函数f加上一些高斯噪声：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file223.jpg

如果f是X的线性函数，那么这个假设本质上与我们在第四章讨论简单线性回归时使用的假设相同。在本章中，我们将通过对f设置先验来使用一个更一般的表达式。通过这种方式，我们将能够得到比线性更复杂的函数。如果我们决定使用高斯过程作为这个先验，那么我们可以写成：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file224.jpg

这里，https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/GP.PNG表示一个均值函数μ[X]和协方差函数K(X,X′)的高斯过程。尽管在实践中，我们总是处理有限对象，但我们使用函数这个词来表示从数学上讲，均值和协方差是无限对象。

我之前提到过，使用高斯分布非常方便。例如，如果先验分布是一个高斯过程，而似然是高斯分布，那么后验分布也是高斯过程，并且我们可以解析地计算它。此外，使用高斯似然的好处是我们可以边际化高斯过程，这大大减少了我们需要从中抽样的参数空间的大小。PyMC 中的高斯过程模块利用了这一点，并为高斯和非高斯似然提供了不同的实现方式。在接下来的部分中，我们将探索这两者。

8.6 使用 PyMC 进行高斯过程回归

图 8.4中的灰色线是一个正弦函数。我们假设我们不知道这个函数，而是只有一组数据点（点）。然后我们使用高斯过程来逼近生成这些数据点的函数。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file225.png

图 8.4：从已知函数（线）生成的合成数据（点）

高斯过程在 PyMC 中实现为一系列 Python 类，这些类与我们之前见过的模型稍有不同；然而，代码依然非常PyMConic。我在接下来的代码中添加了一些注释，帮助你理解在 PyMC 中定义高斯过程的关键步骤。

代码 8.3

# A one-dimensional column vector of inputs. 
X = x[:, None] 

with pm.Model() as model_reg: 
    # hyperprior for lengthscale kernel parameter 
    ℓ = pm.InverseGamma("ℓ", 7, 17) 
    # instanciate a covariance function 
    cov = pm.gp.cov.ExpQuad(1, ls=ℓ) 
    # instanciate a GP prior 
    gp = pm.gp.Marginal(cov_func=cov) 
    σ = pm.HalfNormal('σ', 25) 
    # Class representing that the observed data is a GP plus Gaussian noise 
    y_pred = gp.marginal_likelihood('y_pred', X=X, y=y, sigma=σ) 

    idata_reg = pm.sample()

请注意，我们没有使用高斯似然，而是使用了gp.marginal_likelihood方法。这个方法利用了后验分布有封闭形式这一事实，正如前面部分所解释的那样。

好的，现在我们已经计算了后验分布，接下来让我们看看如何得到均值拟合函数的预测。我们可以通过使用gp.conditional计算在新输入位置上评估的条件分布来实现这一点。

代码 8.4

X_new = np.linspace(np.floor(x.min()), np.ceil(x.max()), 100)[:,None] 
with model_reg: 
    f_pred = gp.conditional('f_pred', X_new)

结果，我们得到了一个新的 PyMC 随机变量f_pred，我们可以用它来从后验预测分布中获取样本：

代码 8.5

with model_reg: 
    idata_subset = idata_reg.sel(draw=slice(0, None, 100)) 
    pred_samples = pm.sample_posterior_predictive(idata_subset, 
                                                  var_names=["f_pred"]) 

f_pred = (pred_samples. 
         posterior_predictive.stack(samples=("chain", "draw"))['f_pred']. 
         values)

现在我们可以在原始数据上绘制拟合函数，以直观检查它们与数据的拟合程度以及我们预测的相关不确定性。正如我们在第四章中对线性模型所做的那样，我们将展示几种不同的方式来绘制相同的结果。图 8.5展示了拟合函数的线条。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file226.png

图 8.5：线条表示从model_reg的后验均值中抽取的样本

或者，我们可以使用辅助函数pm.gp.util.plot_gp_dist来获得如图 8.6所示的漂亮图表。在这个图表中，每条带状线表示不同的百分位数，从 99 百分位（浅灰色）到 51 百分位（深灰色）。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file227.png

图 8.6：使用plot_gp_dist函数绘制的从model_reg的后验样本

另一种选择是计算在参数空间中给定点处的条件分布的均值向量和标准差。在图 8.7中，我们使用均值（基于轨迹中的样本）来表示ℓ和 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/e.png。我们可以使用gp.predict方法来计算均值和方差。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file228.png

图 8.7：model_reg的后验均值，带有 1 和 2 标准差的带宽

正如我们在第 4章中看到的，我们可以使用一个具有非高斯似然度和适当逆链接函数的线性模型来扩展有用线性模型的范围。对于高斯过程（GP），我们也可以做同样的事情。例如，我们可以使用一个带有指数逆链接函数的泊松似然度。对于这样的模型，后验分布不再是解析可解的，但我们仍然可以使用数值方法来近似计算。在接下来的部分中，我们将讨论这些类型的模型。

8.6.1 设置长度尺度的先验

对于长度尺度参数，避免为零的先验通常效果更好。正如我们之前看到的，ℓ控制着函数的平滑度，因此，ℓ为 0 意味着函数不平滑；我们会得到类似于图 8.1中的“第一个”函数。但更重要的原因是，对于大于 0 但仍低于协变量最小间距的ℓ值，可能会出现一些不良效应。本质上，在该点以下，似然度无法区分不同的长度尺度，因此它们都同样好。这是一个不可识别性问题。因此，我们将得到一个倾向于过拟合并且准确插值输入数据的高斯过程（GP）。此外，MCMC 采样器会遇到更多困难，可能会出现更长的采样时间或简单的不可靠样本。类似的情况发生在数据范围之外。如果数据范围为 10 且ℓ> = 10，这意味着函数是平坦的。再往后，您（以及似然度）无法区分不同的参数值。因此，即使您不确定函数是多么平滑或起伏，您仍然可以设置一个先验，避免极低或极高的ℓ值。例如，为了获得先验pm.InverseGamma("ℓ", 7, 17)，我们请求 PreliZ 提供一个最大熵先验，其中 0.95 的质量分布在 1 到 5 之间：

代码 8.6

pz.maxent(pz.InverseGamma(), 1, 5, 0.95)

逆伽马分布（InverseGamma）是一个常见的选择。与伽马分布（Gamma）类似，它允许我们设置一个避免为 0 的先验，但与伽马分布不同的是，逆伽马分布在接近 0 时有一个较轻的尾部，换句话说，它为小值分配的质量较少。

本章剩余部分，我们将使用函数get_ig_params从协变量的尺度中获取弱信息先验。你可以在附带的代码中找到详细信息（github.com/aloctavodia/BAP3），但本质上，我们使用 PreliZ 中的maxent函数，将大部分先验质量设置在一个与协变量范围兼容的区间内。

8.7 高斯过程分类

在第四章中，我们看到了如何使用线性模型来分类数据。我们使用了伯努利似然和逻辑逆链接函数。然后，我们应用了边界决策规则。在这一节中，我们将做同样的事情，但这次使用 GP 而不是线性模型。就像我们在第四章中的model_lrs一样，我们将使用包含两类setosa和versicolor、一个预测变量sepal length的鸢尾花数据集。

对于这个模型，我们不能使用pm.gp.Marginal类，因为该类仅限于高斯似然，因为它利用了 GP 先验与高斯似然组合的数学可处理性。相反，我们需要使用更通用的pm.gp.Latent类。

代码 8.7

with pm.Model() as model_iris: 
    ℓ = pm.InverseGamma('ℓ', *get_ig_params(x_1)) 
    cov = pm.gp.cov.ExpQuad(1, ℓ) 
    gp = pm.gp.Latent(cov_func=cov) 
    f = gp.prior("f", X=X_1) 
    # logistic inverse link function and Bernoulli likelihood 
    y_ = pm.Bernoulli("y", p=pm.math.sigmoid(f), observed=y) 
    idata_iris = pm.sample()

如我们所见，图 8.8看起来与图 4.11非常相似。请花些时间对比这两张图。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file229.png

图 8.8：逻辑回归，model_lrs的结果

你可能已经注意到，推断出的函数看起来类似于 sigmoid 曲线，除了在较低的sepal_length值时尾部向上，较高的sepal_length值时尾部向下。为什么会这样呢？因为当数据很少或没有数据时，GP 后验倾向于恢复到 GP 先验。如果我们考虑到在没有数据的情况下，后验基本上变成了先验，这就很有意义。

如果我们唯一关注的是决策边界，那么尾部的行为可能是无关紧要的。但如果我们想要建模在不同sepal_length值下属于 setosa 或 versicolor 的概率，我们应该做一些事情来改善尾部的模型。实现这一目标的一种方法是为高斯过程添加更多结构。高斯过程的一个非常好的特点是，我们可以组合协方差函数。因此，对于下一个模型，我们将组合三个核函数：指数二次核、线性核和白噪声核。

线性核的效果是在数据边界处使尾部趋近于 0 或 1。此外，我们使用白噪声核作为稳定协方差矩阵计算的一种技巧。高斯过程的核函数受到限制，以保证生成的协方差矩阵是正定的。然而，数值误差可能导致违反这一条件。这种问题的表现之一是，在计算拟合函数的后验预测样本时会得到 NaN。缓解此错误的一种方法是通过添加噪声来稳定计算。事实上，PyMC 在后台已经做了类似的处理，但有时需要更多的噪声，如下面的代码所示：

代码 8.8

with pm.Model() as model_iris2: 
    ℓ = pm.InverseGamma('ℓ', *get_ig_params(x_1)) 
    c = pm.Normal('c', x_1.min()) 
    τ = pm.HalfNormal('τ', 5) 
    cov = (pm.gp.cov.ExpQuad(1, ℓ) + 
           τ * pm.gp.cov.Linear(1, c) + 
           pm.gp.cov.WhiteNoise(1E-5)) 
    gp = pm.gp.Latent(cov_func=cov) 
    f = gp.prior("f", X=X_1) 
    # logistic inverse link function and Bernoulli likelihood 
    y_ = pm.Bernoulli("y", p=pm.math.sigmoid(f), observed=y) 
    idata_iris2 = pm.sample()

我们可以在图 8.9中看到这个模型的结果。注意，这个图看起来现在与图 4.11相比，更加相似，而不是图 8.8。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file230.png

图 8.9：逻辑回归，model_lrs的结果

本节讨论的示例有两个主要目的：

• 显示如何轻松地将核函数组合在一起，以获得更具表现力的模型

• 显示如何使用高斯过程恢复逻辑回归

关于第二点，逻辑回归确实是高斯过程的一个特例，因为简单的线性回归只是高斯过程的一个特例。事实上，许多已知的模型可以看作是高斯过程的特例，或者至少它们与高斯过程有某种联系。如果你想了解更多，可以阅读 Kevin Murphy 的《机器学习：一种概率视角》第一版第十五章[Murphy，2012]，以及第二版的第十八章[Murphy，2023]。

8.7.1 高斯过程用于空间流感

实际上，使用高斯过程来建模我们可以通过逻辑回归解决的问题并没有太大意义。相反，我们希望使用高斯过程来建模那些无法通过不太灵活的模型很好捕捉的更复杂的数据。例如，假设我们想要将患病的概率作为年龄的函数进行建模。结果表明，年轻人和年老的人比中年人有更高的风险。数据集space_flu.csv是一个受上述描述启发的合成数据集。图 8.10展示了它的图形。

让我们拟合以下模型并绘制结果：

代码 8.9

with pm.Model() as model_space_flu: 
    ℓ = pm.InverseGamma('ℓ', *get_ig_params(age)) 
    cov = pm.gp.cov.ExpQuad(1, ℓ) + pm.gp.cov.WhiteNoise(1E-5) 
    gp = pm.gp.Latent(cov_func=cov) 
    f = gp.prior('f', X=age) 
    y_ = pm.Bernoulli('y', p=pm.math.sigmoid(f), observed=space_flu) 
    idata_space_flu = pm.sample()

请注意，如图 8.10所示，高斯过程可以很好地拟合这个空间流感数据集，即使数据要求函数比逻辑回归更复杂。对于简单的逻辑回归来说，拟合这些数据几乎是不可能的，除非我们做一些特殊的修改来帮助它（有关这种修改的讨论请参见本章末的练习 6）。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file231.png

图 8.10：逻辑回归，model_space_flu的结果

8.8 Cox 过程

现在我们要建模计数数据。我们将看到两个例子；一个是具有时间变化率的，另一个是具有二维空间变化率的。为此，我们将使用泊松似然，并通过高斯过程来建模该比率。由于泊松分布的比率限制为正值，我们将使用指数作为逆链接函数，正如我们在第四章中的负二项回归中所做的那样。

我们可以将泊松过程看作是在给定空间中点集的分布，其中每个有限的点集都是泊松分布。当泊松过程的比率本身是一个随机过程时，例如一个高斯过程，那么我们就有了 Cox 过程。

8.8.1 煤矿灾难

第一个例子被称为煤矿灾难。这个例子记录了英国从 1851 年到 1962 年的煤矿灾难数据。灾难的数量被认为受到了这一时期安全法规变化的影响。我们希望将灾难的发生率建模为时间的函数。我们的数据集包含一列数据，每条数据对应一次灾难发生的时间。我们将用于拟合数据的模型形式如下：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/Formula_01.PNG

如您所见，这是一个泊松回归问题。此时，您可能会想，我们如何仅凭一个只有灾难日期的单列数据来进行回归分析。答案是将数据离散化，就像我们构建直方图一样。我们将使用箱子的中心作为X变量，箱子内的计数作为Y变量：

代码 8.10

# discretize data 
years = int((coal_df.max() - coal_df.min()).iloc[0]) 
bins = years // 4 
hist, x_edges = np.histogram(coal_df, 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 disasters per year 
y_data = hist

现在我们使用 PyMC 来定义并求解模型：

代码 8.11

with pm.Model() as model_coal: 
    ℓ = pm.InverseGamma('ℓ', *get_ig_params(x_edges)) 
    cov = pm.gp.cov.ExpQuad(1, ls=ℓ) + pm.gp.cov.WhiteNoise(1E-5) 
    gp = pm.gp.Latent(cov_func=cov) 
    f = gp.prior('f', X=x_data) 
    y_pred = pm.Poisson('y_pred', mu=pm.math.exp(f), observed=y_data) 
    idata_coal = pm.sample()

图 8.11 显示了灾难率随时间变化的中位数（白线）。带状区域描述了 50%的 HDI（较暗）和 94%的 HDI（较亮）。底部的黑色标记表示每次灾难发生的时刻。正如我们所看到的，事故率随着时间的推移而下降，除了最初的短暂上升。PyMC 文档包括了煤矿灾难的案例，但从不同的角度进行建模。我强烈建议您查看这个例子，因为它本身非常有用，并且与我们刚刚实现的方法做对比也很有价值。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file232.png

图 8.11：逻辑回归，model_coal的结果

请注意，即使我们对数据进行了分箱，结果仍然是一个平滑的曲线。从这个角度来看，我们可以将 model_coal（以及一般而言，这种类型的模型）视为构建直方图后进行平滑处理。

8.8.2 红木

让我们将刚才使用的方法应用于一个二维空间问题。我们将使用如图 8.12所示的红木数据集。该数据集（采用 GPL 许可证分发）来自 GPstuff 包。数据集包含了某一地区红木树的位置。推断的动机是获得一个速率图，表示某一地区内树木的数量。

与煤矿灾难示例类似，我们需要对数据进行离散化处理：

代码 8.12

# discretize spatial data 
bins = 20 
hist, x1_edges, x2_edges = np.histogram2d( 
    rw_df[1].values, rw_df[0].values, bins=bins) 
# compute the location of the centers of the discretized data 
x1_centers = x1_edges[:-1] + (x1_edges[1] - x1_edges[0]) / 2 
x2_centers = x2_edges[:-1] + (x2_edges[1] - x2_edges[0]) / 2 
# arrange xdata into proper shape for GP 
x_data = [x1_centers[:, None], x2_centers[:, None]] 
# arrange ydata into proper shape for GP 
y_data = hist.flatten().astype(int)

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file233.png

图 8.12：红木数据

请注意，与其做网格化处理，我们将x1和x2数据视为独立的数组。这使我们能够为每个坐标独立地构建协方差矩阵，从而有效地减少计算高斯过程所需的矩阵大小。然后我们使用LatentKron类将两个矩阵结合起来。需要强调的是，这不是一种数值技巧，而是此类矩阵结构的数学特性，因此我们并没有在模型中引入任何近似或误差。我们只是以一种方式表达它，从而实现更快的计算：

代码 8.13

with pm.Model() as model_rw: 
    ℓ = pm.InverseGamma('ℓ', *get_ig_params(x_data), shape=2) 
    cov_func1 = pm.gp.cov.ExpQuad(1, ls=ℓ[0]) 
    cov_func2 = pm.gp.cov.ExpQuad(1, ls=ℓ[1]) 

    gp = pm.gp.LatentKron(cov_funcs=[cov_func1, cov_func2]) 
    f = gp.prior('f', Xs=x_data) 

    y = pm.Poisson('y', mu=pm.math.exp(f), observed=y_data) 
    idata_rw = pm.sample()

在图 8.13中，灰色的阴影越深，树木的生长速率越高。我们可以想象，我们的目标是找到高生长速率的区域，因为我们可能对木材从火灾中恢复的情况感兴趣，或者我们可能对土壤的一些特性感兴趣，并用树木作为代理。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file234.png

图 8.13：逻辑回归，model_rw的结果

8.9 带有空间自相关的回归分析

以下示例摘自Statistical Rethinking: A Bayesian Course with Examples in R and STAN, Second Edition by Richard McElreath, Copyright (2020) by Chapman and Hall/CRC. 经 Taylor & Francis Group 授权复制。我强烈推荐阅读这本书，因为你会发现很多像这样的好例子以及非常好的解释。唯一的警告是，这些书中的例子是用 R/Stan 实现的，但别担心并继续采样；你可以在github.com/pymc-devs/pymc-resources资源中找到这些例子的 Python/PyMC 版本。

对于这个示例，我们有 10 个不同的岛屿社会；对于每一个社会，我们都有它们使用的工具数量。一些理论预测，较大的人口比小人口能发展和维持更多的工具。因此，我们有一个回归问题，其中因变量是工具的数量，独立变量是人口。因为工具数量是计数变量，所以我们可以使用泊松分布。此外，我们有充分的理论依据认为，人口的对数比绝对人口大小更合适，因为真正重要的（根据理论）是人口的数量级。

到目前为止，我们所考虑的模型是泊松回归，但这里有个有趣的部分。另一个影响工具数量的重要因素是岛屿社会之间的接触率。将接触率纳入我们模型的一种方法是收集这些社会在历史上接触的频率信息，并创建一个分类变量，如低/高接触率。另一种方法是使用社会之间的距离作为接触率的代理，因为合理的假设是，地理上较近的社会比远离的社会更频繁接触。

工具数量、人口规模和坐标存储在本书 GitHub 仓库中的islands.csv文件里（github.com/aloctavodia/BAP3）。

忽略先验分布，我们将要构建的模型是：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/Formula_02.PNG

这个模型是一个线性模型加上一个高斯过程（GP）项。我们用线性部分来建模人口对数的影响，用高斯过程项来建模距离/接触率的影响。通过这种方式，我们实际上在有效地纳入技术暴露的相似性度量（从距离矩阵估算）。因此，我们不会假设总数量仅仅是人口的结果，且各社会间相互独立，而是将每个社会的工具数量建模为其空间分布的函数。

关于空间分布的信息是以纬度和经度的形式存在的，但 PyMC 中的核函数假设所有距离都是欧几里得距离。这可能会带来问题。绕过这个问题的最简洁方法可能是使用考虑岛屿处于大致球形地球上的距离。例如，我们可以使用哈弗辛距离，它基于经纬度计算两个点之间的大圆距离。大圆距离是球面上两点之间的最短距离，是沿着球面测量的。为了使用这种距离，我们需要创建一个新的核函数，如下一个代码块所示。如果你对 Python 中的类不太熟悉，你只需要知道我所做的是复制了 PyMC 代码库中的ExpQuad代码，并稍微修改它，创建了一个新的类ExpQuadHaversine。最大变化是添加了函数/方法haversine_distance。

代码 8.14

class ExpQuadHaversine(pm.gp.cov.Stationary): 
    def __init__(self, input_dims, ls, ls_inv=None, r=6371, active_dims=None): 
        super().__init__(input_dims, ls=ls, ls_inv=ls_inv, active_dims=active_dims) 
        self.r = r # earth radius in km 

    def haversine_distance(self, X): 
        lat = np.radians(X[:, 0]) 
        lon = np.radians(X[:, 1]) 
        latd = lat[:,None] - lat 
        lond = lon[:,None] - lon 
        d = pt.cos(lat[:,None]) * pt.cos(lat) 
        a = pt.sin(latd / 2)** 2 + d * pt.sin(lond / 2)** 2 
        c = 2 * pt.arctan2(pt.sqrt(a), pt.sqrt(1 - a)) 
        return self.r * c 

    def full(self, X, _): 
        return pt.exp(-0.5 * self.haversine_distance(X)**2)

现在我们已经定义了类ExpQuadHaversine，可以像使用之前的内置核函数那样，使用它来定义协方差矩阵。对于这个模型，我们将引入另一个变化。我们将定义一个参数η。这个参数的作用是对高斯过程在 y 轴方向进行缩放。通常，我们会为高斯过程定义ℓ和η这两个参数。

代码 8.15

with pm.Model() as model_islands: 
    η = pm.Exponential('η', 2) 
    ℓ = pm.InverseGamma('ℓ', *get_ig_params(islands_dist)) 

    cov = η * ExpQuadHaversine(2, ls=ℓ) 
    gp = pm.gp.Latent(cov_func=cov) 
    f = gp.prior('f', X=X) 

    *α* = pm.Normal('*α*', 0, 5) 
    *β* = pm.Normal('*β*', 0, 1) 
    μ = pm.math.exp(*α* + *β* * log_pop + f) 
    _ = pm.Poisson('tt_pred', μ, observed=total_tools) 

    idata_islands = pm.sample()

为了理解协方差函数相对于距离的后验分布，我们可以绘制一些来自后验分布的样本，如 图 8.14 所示。黑色曲线表示每个距离的后验中位数协方差，灰色曲线表示从 ℓ 和 η 的联合后验分布中采样的函数。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file235.png

图 8.14：空间协方差的后验分布

图 8.14 中的粗黑线表示社会对之间的协方差的后验中位数与距离的关系。我们使用中位数是因为 ℓ 和 η 的分布非常偏斜。我们可以看到，协方差的平均值并不高，而且在大约 2,000 公里时几乎降到 0。细线表示不确定性，我们可以看到存在很大的不确定性。

现在让我们看看根据模型和数据，岛屿社会之间的相关性有多强。为此，我们必须将协方差矩阵转换为相关矩阵。详情请参见附带的代码。图 8.15 显示了均值相关矩阵的热图。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file236.png

图 8.15：后验均值相关矩阵

从其他观察中脱颖而出的两个观点是，首先，夏威夷非常孤独。这是有道理的，因为夏威夷距离其他岛屿社会非常遥远。其次，我们可以看到马莱库拉（Ml）、提科皮亚（Ti）和圣克鲁兹（SC）之间高度相关。这也是合理的，因为这些社会彼此非常接近，并且它们也有相似数量的工具。

图 8.16 的左面板本质上是一个地图。岛屿社会被表示为它们相对的位置。线条表示社会之间的后验中位数相关性。线条的透明度与相关性的值成正比。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file237.png

图 8.16：空间协方差的后验分布

在图 8.16的右侧面板中，我们再次展示了后验中位数相关性，但这次是根据对数人口与工具总数的关系来绘制的。虚线表示工具的中位数和 94%的高密度区间（HDI）作为对数人口的函数。在图 8.16的两个面板中，点的大小与每个岛屿社会的人口成正比。注意，马列库拉、提科皮亚和圣克鲁斯之间的相关性描述了它们工具数量相对较少，接近中位数或低于根据其人口预期的工具数量。类似的情况发生在特罗布里安群岛和马努斯岛；它们地理位置接近，且工具数量低于预期。汤加的工具数量远高于其人口预期，并且与斐济有较高的相关性。从某种程度上讲，模型告诉我们，汤加对斐济的影响是积极的，增加了工具的总数，同时抵消了对其邻近地区马列库拉、提科皮亚和圣克鲁斯的影响。

8.10 希尔伯特空间高斯过程（HSGP）

高斯过程可能会比较慢。主要原因是它们的计算需要我们对一个矩阵进行求逆，而该矩阵的大小会随着观察数据的增多而增长。这个操作计算成本较高，而且扩展性不好。因此，围绕高斯过程的研究大部分集中在寻找近似方法，以便更快地计算它们并使其能够扩展到大规模数据。

我们将只讨论其中一种近似方法，即希尔伯特空间高斯过程（HSGP），而不深入探讨这种近似是如何实现的。从概念上讲，我们可以将其视为一种基函数展开，类似于样条函数的构建方式（参见第六章）。这种近似的结果是，它将矩阵求逆操作转化为矩阵乘法，这是一种速度更快的操作。

但什么时候它能起作用？

我们只能在低维度（大约 1 到 3 或 4 维）中使用 HSGP，并且仅限于某些核函数，如指数二次核或马特恩核。原因是，为了使 HSGP 近似法生效，核函数必须以一种特殊的形式表示，即功率谱密度形式，而并非所有核函数都可以用这种形式表示。

在 PyMC 中使用 HSGP 近似方法是非常直接的，我们将在自行车数据集上演示这一点。我们希望将出租自行车的数量建模为一天中小时数的函数。以下代码块展示了此模型的 PyMC 实现。

代码 8.16

with pm.Model() as model_hsgp: 
    ℓ = pm.InverseGamma('ℓ', *get_ig_params(X)) 

    cov = pm.gp.cov.ExpQuad(1, ls=ℓ) 
    gp = pm.gp.HSGP(m=[10], c=1.5, cov_func=cov) 

    f = gp.prior('f', X=X) 
    *α* = pm.HalfNormal('*α*', 1) 
    _ = pm.NegativeBinomial("obs", np.exp(f), *α*, observed=y) 

    idata_hsgp = pm.sample()

与之前的高斯过程模型的主要区别在于，使用了pm.gp.HSGP(.)类，而不是用于非高斯似然和标准高斯过程的pm.gp.Latent(.)类。pm.gp.HSGP(.)类有两个参数：

• m 是我们用来逼近 GP 的基本函数数量。m 值越大，逼近效果越好，但计算成本也越高。

• c 是一个边界因子。对于一个固定且足够大的 m 值，c 主要影响均值函数在边界附近的逼近。它不应小于 1.2（如果使用较小的值，PyMC 会给出警告），通常 1.5 是一个不错的选择。改变此参数不会影响计算速度。

我们设置 m=10，部分原因是我们喜欢十进制系统，部分原因是基于 Riutort-Mayol et al. 发表的论文《Practical Hilbert space approximate Bayesian Gaussian processes for probabilistic programming》中的建议 [2022]。实际上，只要 m 和 c 的值在某个范围内，与长度尺度的先验一致，结果对它们的精确值是稳健的。关于 HSGP 如何工作以及如何在实践中使用它的一些建议，可以参考 Riutort-Mayol et al. [2022]。

现在让我们看看结果。图 8.17 显示了黑色的均值后验 GP 和来自 GP 后验的 100 个样本（灰色线）。你可以将这些结果与使用样条得到的结果进行比较（见图 6.8）。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file238.png

图 8.17：HSGP 模型对租赁自行车的后验均值，作为一天中时间的函数

HSGP 近似也已在 Bambi 中实现。让我们看看如何使用它。

8.10.1 HSGP 与 Bambi

要使用 Bambi 拟合之前的模型，我们需要编写以下代码：

代码 8.17

bmb.Model("rented ∼ 0 + hsgp(hour, m=10, c=1.5)", bikes, 
          family="negativebinomial")

这样是可行的，但我们将为 Bambi 提供先验，就像我们在 PyMC 中做的那样。这将导致更快的采样和更可靠的样本。

正如我们在第 6 章中看到的那样，要在 Bambi 中定义先验，我们只需将字典传递给 bmb.Model 的 priors 参数。但我们必须注意，HSGP 项不会接收先验。相反，我们需要为 ℓ（在 Bambi 中称为 ell）和 η（在 Bambi 中称为 sigma）定义先验，并将这些先验传递给 HSGP 项。还有一件事：和之前的模型一样，我们没有使用 η，但由于 Bambi 需要它，我们使用了一个小技巧来定义一个基本上是 1 的先验。

代码 8.18

prior_gp = { 
    "sigma": bmb.Prior("Gamma", mu=1, sigma=0.01), 
    "ell": bmb.Prior("InverseGamma", **get_ig_params(X)) 
} 
priors = { 
    "hsgp(hour, m=10, c=1.5)": prior_gp, 
    "alpha": bmb.Prior("HalfNormal", sigma=1) 
} 

model_hsb = bmb.Model("rented ∼ 0 + hsgp(hour, m=10, c=1.5)", bikes, 
                      family="negativebinomial", 
                      priors=priors) 

idata_hsb = model_hsb.fit()

我邀请你检查 Bambi 计算的参数，它们与我们通过 PyMC 得到的非常相似。图 8.18 显示了黑色的均值后验 GP 和 94% HDI 的带状区域。该图是通过 bmb.interpret.plot_predictions 生成的。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file239.png

图 8.18：HSGP 模型对租赁自行车的后验均值，作为一天中时间的函数，使用 Bambi

在本节中，我们探讨了 HSGP 的概念，作为一种强大的近似方法，可以将高斯过程扩展到大数据集。通过将 PyMC 和 Bambi 的灵活性与 HSGP 提供的可扩展性相结合，研究人员和实践者可以更有效地解决复杂的建模任务，为在日益庞大和复杂的数据集上应用高斯过程铺平道路。

8.11 总结

高斯过程是多元高斯分布的一种推广，适用于无限多个维度，完全由均值函数和协方差函数指定。由于我们可以概念上将函数视为无限长的向量，因此可以将高斯过程作为函数的先验。在实践中，我们并不处理无限对象，而是处理与数据点数量相等维度的多元高斯分布。为了定义其相应的协方差函数，我们使用了适当参数化的核；通过学习这些超参数，我们最终学会了关于任意复杂函数的知识。

在本章中，我们简要介绍了高斯过程（GP）。我们涵盖了回归、半参数模型（岛屿示例）、将两个或多个核组合以更好地描述未知函数，以及高斯过程如何用于分类任务。还有许多其他主题我们可以讨论。然而，我希望这段高斯过程的介绍足以激励你继续使用、阅读和学习高斯过程及贝叶斯非参数模型。

8.12 练习

1. 对于协方差函数和核部分中的示例，确保理解输入数据与生成的协方差矩阵之间的关系。尝试使用其他输入，如data = np.random.normal(size=4)。

2. 重新运行生成图 8.3的代码，并将从高斯过程（GP）先验中获得的样本数量增加到约 200。原始图中的样本数量为 2。生成值的范围是什么？

3. 对于前一个练习中生成的图表，计算每个点的标准差。按照以下形式进行操作：

   • 从视觉上看，仅仅观察图表

   • 直接来自于pz.MVNormal(.).rvs生成的值

   • 通过检查协方差矩阵（如果有疑问，请回到练习 1）

   从这三种方法中得到的值是否匹配？

4. 使用测试点np.linspace(np.floor(x.min()), 20, 100)[:,None]并重新运行model_reg。绘制结果。你观察到了什么？这与 GP 先验的规范有何关系？

5. 重复练习 1，但这次使用线性核（参见附带的线性核代码）。

6. 请查看 PyMC 文档中的www.pymc.io/projects/examples/en/latest/gaussian_processes/GP-MeansAndCovs.html。

7. 对space_flu数据运行逻辑回归模型。你看到了什么？你能解释结果吗？

8. 修改逻辑回归模型以适应数据。提示：使用二阶多项式。

9. 将煤矿灾难的模型与 PyMC 文档中的模型进行比较（www.pymc.io/projects/docs/en/stable/learn/core_notebooks/pymc_overview.html#case-study-2-coal-mining-disasters）。描述这两个模型在模型规范和结果方面的差异。

加入我们的社区 Discord 空间

加入我们的 Discord 社区，结识志同道合的人，并与 5000 多名成员一起学习：packt.link/bayesian

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file1.png

第九章

贝叶斯加法回归树

个人而言，我们是一个水滴。一起时，我们是海洋。—— 佐藤隆之

在上一章中，我们讨论了高斯过程（GPs），一种用于回归的非参数模型。在本章中，我们将学习另一种非参数回归模型，称为贝叶斯加法回归树，或者亲切地称为 BART。我们可以从多个不同的角度来看待 BART。它可以看作是决策树的一个集成，每棵树在整体数据理解中扮演着独特的角色和贡献。这些树在贝叶斯先验的指导下和谐工作，以捕捉数据的细微差别，避免个体过拟合的陷阱。通常，BART 作为一个独立的模型进行讨论，实施该模型的软件通常仅限于一个或少数几个模型。在本章中，我们将采用不同的方式，使用 PyMC-BART，这是一个允许在 PyMC 中使用 BART 模型的 Python 库。

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

• 决策树

• BART 模型

• 使用 BART 的灵活回归

• 部分依赖图

• 单个条件期望图

• 变量选择

9.1 决策树

在深入讨论 BART 模型之前，先花点时间了解一下什么是决策树。决策树就像一个流程图，指导你通过不同的问题，直到你做出最终选择。例如，假设你每天早上需要决定穿什么鞋子。为了做出决定，你可能会问自己一系列问题：“今天暖和吗？”如果是的话，你可能会问更具体的问题，比如“我需要去办公室吗？”最终，你会停止提问并得到一个输出值，比如拖鞋、运动鞋、靴子、莫卡辛鞋等。

该流程图可以方便地编码为树状结构，在树的根部放置更一般性的问题，然后沿着树结构向更具体的问题推进，最终到达树的叶子节点，输出不同类型鞋子的结果。树是计算机科学和数据分析中非常常见的数据结构。

更正式地说，树是节点和连接这些节点的顶点的集合。带有问题的节点称为决策节点，带有树的输出（如鞋子）的节点称为叶子节点。当答案是“是”或“否”时，我们就有了一个二叉树，因为每个节点最多可以有两个子节点。图 9.1 显示了一个决策树。圆角方块是叶子节点，普通方块是决策节点。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file240.png

图 9.1：选择鞋类的决策树。

我们可以使用决策树进行分类，即返回离散类别，如运动鞋、拖鞋、凉鞋等。但我们也可以用它们进行回归，即返回连续结果，如 4.24 或 20.9（以及介于两者之间的任何值）。通常，这些树被称为回归树。图 9.2展示了左侧的回归树。我们还可以将回归树视为分段阶梯函数的表示，如图 9.2右侧所示。这与表示平滑函数（至少在某种程度上）的三次样条或高斯过程（GPs）不同。树可以足够灵活，以提供良好的平滑函数的实用近似。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file241.png

图 9.2：左侧为回归树，右侧为对应的分段阶梯函数

树可以非常灵活；在极端情况下，我们可以有一棵树，其中的叶节点数量与观测值一样多，这棵树将完美拟合数据。正如我们在第五章中看到的，除非我们添加一些正则化，否则这可能不是一个好主意。在贝叶斯术语中，我们可以通过先验来实现这种正则化。例如，我们可以设置一个先验，使树较浅。通过这种方式，我们使得树的节点数与数据点数相等的情况非常不太可能发生。

9.2 BART 模型

贝叶斯加性回归树（BART）模型是m棵树的和，我们用它来逼近一个函数[Chipman et al.，2010]。为了完成模型，我们需要设置树的先验。这些先验的主要作用是防止过拟合，同时保留树所提供的灵活性。先验设计旨在使得每棵树相对较浅，并且叶节点的值相对较小。

PyMC 不直接支持 BART 模型，但我们可以使用 PyMC-BART，一个扩展 PyMC 功能以支持 BART 模型的 Python 模块。PyMC-BART 提供：

• 一个与 PyMC 中其他分布（如pm.Normal、pm.Poisson等）非常相似的 BART 随机变量。

• 一个名为 PGBART 的采样器，因为树不能使用 PyMC 的默认步进方法（如 NUTS 或 Metropolis）进行采样。

• 以下是一些帮助处理 BART 模型结果的实用函数：

  • pmb.plot_pdp：用于生成部分依赖图的函数[Friedman，2001]。

  • pmb.plot_ice：用于生成个体条件期望图的函数[Goldstein et al.，2013]。

  • pmb.plot_variable_importance：用于估算变量重要性的函数。

  • pmb.plot_convergence：一个绘制 BART 随机变量的有效样本大小的经验累积分布以及https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/hat_R.png值的函数。

BART 是步进函数的先验

我们可以将 BART 看作是分段常数函数的先验。此外，当树的数量m → ∞时，BART 会收敛为一个处处不可微的高斯过程。

在接下来的部分，我们将重点讨论 BART 的应用部分，特别是如何使用 PyMC-BART。如果您有兴趣了解更多有关 BART 模型工作原理的细节，PyMC-BART 的实现细节，以及更改 PyMC-BART 的超参数如何影响结果，建议阅读 Quiroga 等人 [2022]。

9.2.1 巴特企鹅

假设出于某种原因，我们有兴趣将企鹅的体重作为其他身体指标的函数进行建模。下面的代码块展示了一个针对这种问题的 BART 模型。在这个示例中，X = "flipper_length", "bill_depth", "bill_length"]，而Y是body_mass：

代码 9.1

with pm.Model() as model_pen: 

    σ = pm.HalfNormal("σ", 1) 

    μ = pmb.BART("μ", X, Y) 

    y = pm.Normal("y", mu=μ, sigma=σ, observed=Y) 

    idata_pen = pm.sample() 

    pm.sample_posterior_predictive(idata_pen, extend_inferencedata=True)

我们可以看到，使用 PyMC-BART 通过 PyMC 定义 BART 模型是很简单的。基本上，我们需要定义一个带有参数X（协变量）和Y（响应变量）的 BART 随机变量。除此之外，模型的其余部分应该与其他回归模型非常相似。像其他回归模型一样，μ的长度将与观察值的数量相同。

从理论上讲，这些树仅仅是X的函数，但 PyMC-BART 需要Y来获取树叶节点上方差初始值的估计。

一旦我们使用 BART 变量拟合了模型，剩余的工作流程就和通常一样。例如，我们可以通过调用az.plot_ppc(.)来计算后验预测检验，结果将类似于图 9.3。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file242.png

图 9.3：model_pen的后验预测检验

图 9.3显示了一个合理的拟合。值得注意的是，即使我们使用正态似然，也没有得到负的质量。然而，使用 PyMC 和 PyMC-BART 时，尝试其他似然分布非常简单；只需将正态分布替换为其他分布（如 Gamma 分布或截断正态分布），就像在常规的 PyMC 模型中那样，就能顺利进行。然后，您可以使用后验预测检验和 LOO，如第五章中讨论的那样。

在接下来的几节中，我们将讨论如何使用和解释 PyMC-BART 提供的工具函数（pmb.plot_convergence除外，该函数将在第十章中与其他诊断方法一起讨论）。

9.2.2 部分依赖图

部分依赖图（PDP）是一种在 BART 文献中广泛使用的图形工具，但并非仅限于 BART。原则上，它可以与任何方法或模型一起使用。它的基本原理是绘制给定协变量 X[i] 下的预测响应，同时对其余协变量 X[−i] 进行平均。因此，实质上，我们是在绘制每个协变量对响应变量的贡献，同时保持其他变量不变。对于 BART 和其他基于树的方法，有一个特点是 PDP 的计算可以在不重新拟合模型到合成数据的情况下进行；相反，它可以通过已拟合的树高效地计算出来。这使得 BART 成为模型可解释性和理解单个特征对预测影响的一个有吸引力的选择。

一旦像 model_pen 这样的模型已经拟合，我们可以通过以下方式计算部分依赖图（PDP）：

代码 9.2

pmb.plot_pdp(μ, X, Y)

请注意，我们传递了 BART 随机变量、协变量和响应变量。响应变量其实并非必需，但如果传递了且它是 pandas Series 类型，它将使用其名称作为 y 轴标签。

图 9.4 显示了来自 model_pen 的部分依赖图示例。我们可以看到，flipper_length 显示出最大的效应，呈线性关系，而另外两个变量则基本保持平坦的响应，表明它们对响应变量的部分贡献不大。对于对响应变量贡献为零的变量，其预期的 PDP 将是一个平坦的常数线，值为响应变量的平均值。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file243.png

图 9.4：model_pen 的部分依赖图

在 图 9.4 中，我们可以看到，最大的贡献来自 flipper_length，但这并不意味着其他两个变量与 body_mass 没有关系。我们只能说，考虑到模型中已经包含了 flipper_length，其他两个变量的影响是最小的。

9.2.3 个体条件图

在计算部分依赖图时，我们假设变量 X[i] 和 X[−i] 之间没有相关性。在许多实际问题中，情况往往并非如此，部分依赖图可能会掩盖数据中的关系。然而，如果所选变量子集之间的依赖关系不太强，则部分依赖图仍然可以作为有用的总结工具 Friedman [2001]。

个体条件期望（ICE）图与部分依赖图（PDP）密切相关。不同之处在于，我们不是绘制目标协变量对预测响应的平均部分效应，而是在给定固定值（默认为 10）下，绘制 n 条条件期望曲线。也就是说，ICE 图中的每条曲线反映了在固定的 X[i] 值下，协变量 X[i] 对预测响应的部分影响。

一旦像 model_pen 这样的模型已经拟合，我们可以通过以下命令计算 ICE 图：

代码 9.3

pmb.plot_ice(μ, X, Y)

这个签名与 PDP 相同。结果显示在图 9.5中。灰色曲线是不同值下的条件期望曲线。如果我们对它们取平均值，我们应该得到 PDP 曲线（黑色）。如果 ICE 图中的曲线大多平行，这意味着协变量对响应变量的贡献大多是独立的。这对于flipper_length和bill_length来说是成立的。在这种情况下，ICE 图和 PDP 图传达的是相同的信息。然而，如果曲线交叉，则表示贡献是非独立的。在这种情况下，PDP 会隐藏这些效应。我们可以在下图中看到bill_depth的例子：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file244.png

图 9.5：model_pen的个体条件期望图

默认情况下，ICE 图是居中的，这意味着灰色曲线围绕在 x 轴最低值处计算的部分响应进行居中。这有助于解释图形：例如，可以更容易地看出线条是否交叉。这也解释了为什么图 9.5中的 y 轴刻度与图 9.4中的刻度不同。你可以通过参数centered=False来修改这一点。

9.2.4 使用 BART 进行变量选择

在第 6章中，我们已经讨论过变量选择，并解释了在某些场景下我们可能有兴趣选择一个变量子集。PyMC-BART 提供了一种非常简单、几乎不需要计算的启发式方法来估计变量重要性。它跟踪每个协变量作为分裂变量被使用的次数。对于 BART 模型，变量重要性是通过在m棵树和所有后验样本上取平均值来计算的。为了进一步简化解释，我们可以报告经过归一化的值，使每个值都位于区间[0, 1]内，且总重要性为 1。

在某些 BART 实现中，变量重要性的估计对树的数量m非常敏感。这些实现的作者建议，在进行变量选择时使用相对较少的树，而在进行模型拟合/预测时使用更多的树。PyMC-BART 不是这种情况，PyMC-BART 的变量重要性估计对树的数量具有鲁棒性。

一旦我们像model_pen这样拟合了一个模型来使用 PyMC-BART 进行变量选择，我们需要做如下操作：

代码 9.4

pmb.plot_variable_importance(idata_pen, μ, X)

请注意，我们传递了推断数据、BART 随机变量和协变量。结果显示在图 9.6中：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file245.png

图 9.6：model_pen的变量重要性图

从顶部面板，我们可以看到flipper_length具有最大的变量重要性值，其次是bill_depth和fill_length。请注意，这与部分依赖图和个体条件期望图在定性上是一致的。

计算一个变量在树中出现的次数这一简单启发式方法存在一些问题。其中一个问题涉及可解释性，因为缺乏一个明确的阈值来区分重要变量和不重要变量，这是一个问题。PyMC-BART 提供了一些帮助。图 9.6的底部面板显示了参考模型生成的预测与使用子模型生成的预测之间的皮尔逊相关系数的平方，参考模型是包含所有协变量的模型，子模型是包含较少协变量的模型。我们可以使用这个图来找到能够做出与参考模型尽可能接近的预测的最小模型。图 9.6告诉我们，只有flipper_length的模型，其预测性能几乎与包含所有三个变量的模型相同。请注意，通过加入bill_depth可能会稍微提高一点性能，但这可能是微不足道的。

现在，让我简要解释一下pmb.plot_variable_importance在背后做了什么。主要有两个近似过程：

• 它不会评估所有可能的协变量组合。相反，它一次添加一个变量，按照变量的重要性顺序进行（图 9.6中的上面子图）。

• 它并不会重新拟合从 1 到 n 个协变量的所有模型。相反，它通过遍历参考模型的后验分布中的树，近似估计去除一个变量的影响，并且剪枝去除没有兴趣变量的分支。这与计算部分依赖图的过程类似，不同之处在于，对于部分依赖图，我们排除了除了一个变量之外的所有变量，而对于变量重要性，我们首先排除除了最重要的一个变量之外的所有变量，然后排除除了最重要的两个变量之外的所有变量，依此类推，直到包括所有变量。

如果这个变量选择过程对你来说听起来很熟悉，那么很可能你一直在关注这一章，并且也看过第六章。这个过程在概念上类似于 Kulprit 所做的工作。在这里，我们也利用了参考模型的概念，并根据其预测分布来评估模型。但相似之处到此为止。PyMC-BART 并没有使用 ELPD，而是使用皮尔逊相关系数的平方，并通过修剪用参考模型拟合的树来估计子模型，而不是通过 Kullback-Liebler 散度投影。

在进入下一个话题之前，让我再补充几点提醒。正如我们在第六章中讨论的那样，使用 Kulprit 的输出时，我们不应过度解读变量的顺序。样本适用于通过pmb.plot_variable_importance生成的图形，如图 9.6。如果两个变量的重要性非常相似，当我们用不同的随机种子重新拟合模型，或数据发生轻微变化（如添加或删除一个数据点）时，顺序可能会发生变化。变量重要性的误差条可以提供帮助，但它们可能低估了真实的变异性。因此，要谨慎对待顺序，将其作为解决问题时的参考。

9.3 分布式 BART 模型

正如我们在第六章中看到的，对于广义线性模型，我们并不限于为均值或位置参数创建线性模型；我们还可以建模其他参数，例如高斯分布的标准差，甚至同时建模均值和标准差。BART 模型也适用同样的原则。

为了举例说明，让我们以自行车数据集为模型。我们将使用rented作为响应变量，hour、temperature、humidity和workday作为预测变量。正如我们之前所做的，我们将使用负二项分布作为似然函数。该分布有两个参数，μ和alpha。我们将为这两个参数使用树的和。以下代码块展示了模型：

代码 9.5

with pm.Model() as model_bb: 
    μ = pmb.BART("μ", X, np.log(Y), shape=(2, 348), separate_trees=True) 
    pm.NegativeBinomial('yl', np.exp(μ[0]), np.exp(μ[1]), observed=Y) 
    idata_bb = pm.sample(2000, 
                         random_seed=123, 
                         pgbart={"batch":(0.05, 0.15)})

让我们花点时间确保理解这个模型。首先，注意我们传递了一个shape参数给pmb.BART()。当separate_trees = True时，这指示 PyMC-BART 拟合两组独立的树和。然后我们索引μ，以便使用第一维度作为负二项分布的μ参数，第二维度作为α参数。如果separate_trees = False，则告诉 PyMC-BART 拟合一组树和，但每棵树在每个叶节点返回 2 个值，而不是 1 个。这样做的好处是算法运行更快，内存使用更少，因为我们只拟合一组树。缺点是模型的灵活性较差。实际上，这两种选项都可以有用，因此你应该选择哪个，取决于你的建模决策。

model_bb的另一个重要方面是我们对μ取指数。这样做是为了确保 NegativeBinomial 分布的μ和α值都为正值。这与我们在广义线性模型中讨论的变换类型相同。PyMC-BART 的特殊之处在于，我们将其逆变换应用于传递给pmb.BART()的Y值。根据我的经验，这有助于 PyMC-BART 找到更好的解。对于具有二项分布或类别分布的模型，无需分别对逻辑回归或 softmax 应用逆变换。PyMC-BART 将二项分布作为特例处理，对于类别分布，我们通过经验发现即使没有逆变换，也能获得良好的结果。需要强调的是，我们传递给pmb.BART()的Y值仅用于初始化 BART 变量的采样。初始化似乎对我们传递的值具有鲁棒性，传递Y或其某种变换在大多数情况下都能很好地工作。

我希望你注意的第三个方面是，我们向pm.sample传递了一个新的参数，即pgbart。这个参数的值是字典"batch":(0.05, 0.15)。为什么要这么做呢？有时，为了获得高质量的样本，必须调整采样器的超参数。在之前的示例中，我们选择省略这一部分，以保持简单和专注。然而，正如我们在第 10章中更深入地讨论的那样，关注这些调整可能变得非常重要。对于 PGBART 采样器的特殊情况，我们可以改变两个超参数。其中一个是num_particles（默认为 10），粒子数量越大，BART 的采样越准确，但代价也越高。另一个是batch；默认值为元组(0.1, 0.1)，意味着在每一步中，采样器在调节阶段拟合 10%的m棵树，在采样阶段也是如此。对于model_bb模型，我们使用了(0.05, 0.15)，即在调节阶段使用 5%（2 棵树），在实际采样阶段使用 15%（7 棵树）。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file246.png

图 9.7：model_bb的部分依赖图

我们可以像图 9.7中一样，探索协变量与响应之间的关系。请注意，变量出现了两次：第一列对应参数μ，第二列对应参数α。我们可以看到，hour对 NegativeBinomial 的两个参数的响应变量影响最大。

9.4 常数和线性响应

默认情况下，PyMC-BART 将拟合在每个叶节点返回单一值的树。这是一种通常非常有效的简单方法。然而，理解其影响是很重要的。例如，这意味着对于任何超出用于拟合模型的观察数据范围的值，预测将是常数。为了验证这一点，回去查看图 9.2。这棵树对于c1以下的任何值都将返回 1.9。请注意，即使我们将多棵树相加，这仍然是事实，因为加总一堆常数值仍然会得到另一个常数值。

这是否是一个问题，取决于你和你应用 BART 模型的上下文。尽管如此，PyMC-BART 提供了一个response参数，可以传递给 BART 随机变量。它的默认值是"constant"。你可以将其更改为"linear"，在这种情况下，PyMC-BART 将在每个叶节点返回一个线性拟合，或者更改为"mix"，这将在采样过程中提出具有常数或线性值的树。

为了举例说明差异，我们来拟合一个非常简单的例子：租赁自行车的数量与温度的关系。观察到的温度值从≈−5 到≈35。拟合此模型后，我们将要求范围为[-20, 45]的样本外后验预测值。因此，我们将设置一个具有可变变量的模型，如第四章所介绍的。

代码 9.6

with pm.Model() as model_tmp1: 
    X_mut1 = pm.MutableData("X_mut1", X) 
    *α* = pm.HalfNormal('*α*', 1) 
    μ = pmb.BART("μ", X_mut1, np.log(Y), m=100, response="linear") 
    _ = pm.NegativeBinomial('yl', np.exp(μ), *α*, observed=Y, shape=μ.shape) 
    idata_tmp1 = pm.sample(random_seed=123)

注意我们将shape=μ.shape传递给了似然函数。这是我们需要做的，以便能够改变X_mut1的形状，这也是 PyMC 的要求，所以你在使用非 BART 模型（如线性回归）时也应该做这件事。

好的，继续这个例子，在随附的代码中，你将找到model_tmp0模型的代码，它与model_tmp1完全相同，只是它具有默认的常数响应。两个模型的结果显示在图 9.8中。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file247.png

图 9.8：常数和线性响应的平均预测

请注意，在数据范围之外（虚线灰色线），具有常数响应的模型的预测确实是常数。哪个模型提供的预测更好？我不确定。我认为，在较低温度范围内，线性响应更好，因为它预测租赁自行车的数量将继续减少，直到最终降至 0。但在较高温度范围内，平台效应或甚至下降应该比上升更可能。我是说，我尝试过在 40 度甚至 42 度时骑自行车，那可不是一个超级愉快的体验。你怎么看？

9.5 选择树的数量

树的数量（m）控制 BART 函数的灵活性。一般来说，默认值 50 应该足以得到一个不错的近似。较大的值，比如 100 或 200，应当提供更精细的结果。通常，增加树的数量不会导致过拟合，因为树的数量越大，叶节点的值越小。

在实际应用中，你可能会担心 m 的值过大，因为 BART 的计算成本（无论是时间还是内存）会随之增加。调节 m 的一种方式是执行 K 折交叉验证，正如 Chipman et al. [2010] 所推荐的那样。另一种选择是通过使用 LOO 来近似交叉验证，正如在第五章中讨论的那样。我们已经观察到 LOO 确实能帮助提供合理的 m 值 [Quiroga et al., 2022]。

9.6 总结

BART 是一种灵活的非参数模型，通过树的总和来近似从数据中获得的未知函数。先验用于规范推理，主要是通过限制树的学习能力，使得没有单一的树能够解释数据，而是树的总和。PyMC-BART 是一个 Python 库，它扩展了 PyMC 以支持 BART 模型。

本章中我们构建了一些 BART 模型，并学习了如何执行变量选择，利用部分依赖图和个体条件图来解释 BART 模型的输出。

9.7 练习

1. 解释以下内容：

   • BART 和线性回归、样条有什么不同？

   • 在什么情况下你可能会选择线性回归而非 BART？

   • 在什么情况下你可能会选择高斯过程而非 BART？

2. 用你自己的话解释为什么多个小树比一棵大树更能拟合模式。两者方法有什么区别？各自有什么权衡？

3. 以下是两个简单的合成数据集。对每个数据集拟合一个 m=50 的 BART 模型。绘制数据和均值拟合函数。描述拟合情况。

   • x = np.linspace(-1, 1., 200) 和 y = np.random.normal(2*x, 0.25)

   • x = np.linspace(-1, 1., 200) 和 y = np.random.normal(x**2, 0.25)

   • 创建你自己的合成数据集。

4. 创建以下数据集 Y = 10sin(πX[0]X[1])+20(X[2] −0*.5)² +10X*[3] +5X[4] + https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/e.png，其中 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/e.png ∼https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/N.PNG(0*,1)，且 X[0:9] ∼https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/U.PNG(0,*1)。这就是所谓的 Friedman 的五维函数。注意，我们实际上有 10 个维度，但最后 5 个是纯噪声。

   • 将 BART 模型拟合到这个数据中。

   • 计算 PDP 和变量重要性（VI）。

   • PDP 和 VI 是否在定性上相符？如何相符？

5. 使用 BART 模型和企鹅数据集。将bill_length、flipper_length、bill_depth、bill_length和body_mass作为协变量，物种Adelie和Chistrap作为响应变量。尝试不同的m值——10、20、50 和 100。使用 LOO 选择合适的值。

6. 查看前一个问题中模型的变量重要性。将结果与使用 Kulprit 构建的同一协变量和响应变量的广义线性模型的结果进行比较，该模型使用 Bambi 构建。

加入我们的社区 Discord 空间

加入我们的 Discord 社区，结识志同道合的人，并与超过 5000 名成员一起学习： packt.link/bayesian

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file1.png

第十章

推断引擎

第一个原则是，你绝不能欺骗自己——而你最容易被自己欺骗。——理查德·费曼

到目前为止，我们主要关注模型构建、结果解释和模型批评。我们已经依赖 pm.sample 函数的魔力为我们计算后验分布。现在我们将集中学习该函数背后推断引擎的一些细节。

概率编程工具的整个目的，比如 PyMC，就是用户不必关心采样是如何进行的，但理解我们如何从后验分布中获取样本对完全理解推断过程非常重要，也有助于我们了解这些方法在何时何种情况下失败，以及如何应对。如果你不关心这些方法是如何工作的，你可以跳过本章的大部分内容，但我强烈建议你至少阅读诊断样本这一节，因为这一节提供了一些指导方针，帮助你检查后验样本是否可靠。计算后验分布的方法有很多。在本章中，我们将讨论一些通用思想，并重点介绍 PyMC 中实现的最重要的方法。我们将学习：

• 推断引擎

• 美特罗波利斯-哈斯廷斯算法

• 哈密尔顿蒙特卡洛

• 序贯蒙特卡洛

• 诊断样本

10.1 推断引擎

虽然贝叶斯方法在概念上简单，但在数学和数值计算上具有挑战性。主要原因是边际似然性，贝叶斯定理中的分母，通常表现为一个难以求解或计算上开销较大的积分。因此，后验通常是通过数值方法使用马尔可夫链****蒙特卡洛（MCMC）家族的算法来估计的。这些方法有时被称为推断引擎，因为至少在理论上，它们能够近似任何概率模型的后验分布。尽管在实际应用中推断并不总是那么理想，但这些方法的存在促使了概率编程语言（如 PyMC）的发展。

概率编程语言的目标是将模型构建过程与推理过程分离，以便于模型构建、评估以及模型修改/扩展的迭代步骤。通过将推理过程（而非模型构建过程）视为黑箱，像 PyMC 这样的概率编程语言用户可以专注于他们的具体问题，而将计算细节交给 PyMC 来处理。这正是我们迄今为止所做的事情。因此，你可能会偏向于认为这是显而易见或自然的方法。但需要注意的是，在概率编程语言出现之前，使用概率模型的人员也习惯于编写自己的采样方法，这些方法通常是根据他们的模型量身定制的，或者他们习惯于简化模型，使其适应某些数学近似。事实上，在某些学术圈子里，这种情况仍然存在。这种量身定制的方法可能更优雅，甚至可以为非常特定的模型提供更高效的后验计算方法，但它也容易出错且耗时，即使对于专家来说也是如此。此外，这种量身定制的方法并不适用于大多数希望利用概率模型解决问题的实践者。像 PyMC 这样的软件邀请来自不同背景的人们使用概率模型，降低了数学和计算的入门门槛。我个人认为这非常棒，也是在邀请我们更多地了解统计建模中的最佳实践，以避免自欺欺人。

之前的章节主要介绍了贝叶斯建模的基础知识；现在我们将从概念层面学习，自动推理是如何实现的，何时以及为什么会失败，失败时该怎么办。

然而，在讨论 MCMC 方法之前，让我先解释另外两种方法，这些方法有时也很有用，并且提供一种直观的理解，说明为什么我们通常使用 MCMC 作为通用方法。

10.2 网格法

网格法是一种简单的暴力法。即使你无法计算整个后验分布，你也许能够逐点计算先验和似然性；这是一个相当常见的情况，甚至可以说是最常见的情形。

假设我们想计算一个单参数模型的后验分布。网格逼近方法如下：

1. 为参数定义一个合理的区间（先验应该会给你一些提示）。

2. 在该区间上放置一个点的网格（通常是等距的）。

3. 对于网格中的每个点，计算似然性与先验的乘积。

可选地，我们可以对计算值进行归一化，即将posterior数组中的每个值除以曲线下的总面积，确保总面积为 1。

以下代码块实现了硬币翻转模型的网格法：

代码 10.1

def posterior_grid(grid_points=50, heads=6, tails=9): 
    """ 
    A grid implementation for the coin-flipping problem 
    """ 
    grid = np.linspace(0, 1, grid_points) 
    prior = np.repeat(1/grid_points, grid_points)  # uniform prior 
    likelihood = pz.Binomial(n=heads+tails, p=grid).pdf(heads) 
    posterior = likelihood * prior 
    posterior /= posterior.sum() * (1/grid_points) 
    return grid, posterior

图 10.1展示了我们在均匀先验下，抛硬币 13 次并观察到 3 次正面的后验。由于我们仅使用了 10 个点的网格，曲线显得非常崎岖。如果增加点的数量，曲线会变得更加平滑，计算结果也会更准确，但成本也会更高。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file248.png

图 10.1：使用网格方法计算的后验

网格方法的最大弊端是，当参数的数量（也称为维度）增多时，这种方法的扩展性较差。我们可以通过一个简单的例子来说明这一点。假设我们想要采样一个单位区间（参见图 10.2），就像抛硬币问题一样，我们使用四个等距的点；这意味着分辨率为 0.25 单位。现在，假设我们有一个二维问题（图 10.2中的正方形），并且我们想用相同分辨率的网格；我们需要 16 个点。最后，对于三维问题，我们需要 64 个点（参见图 10.2中的立方体）。在这个例子中，我们需要 16 倍的资源来从一个边长为 1 的立方体中采样，相比于从一个边长为 1、分辨率为 0.25 的线段中采样。如果我们决定改为分辨率为 0.1 单位，我们将需要为线段采样 10 个点，为立方体采样 1,000 个点。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file249.png

图 10.2：1 维、2 维和 3 维中具有相同分辨率的网格

除了点的数量增长速度外，还有一个现象，并非网格方法或任何其他方法的特性，而是高维空间的特性。随着参数数量的增加，后验概率大部分集中在参数空间的某个区域，而这个区域相对于采样体积来说变得越来越小。这是一个普遍存在的现象，通常被称为“维度灾难”，或者数学家们更喜欢称之为“测度集中”。

维度灾难是指一些相关现象，这些现象在低维空间中不存在，但在高维空间中存在。以下是一些这些现象的例子：

• 随着维度的增加，任何两个样本之间的欧几里得距离趋向于类似其他样本对之间的距离。也就是说，在高维空间中，大多数点之间的距离基本上是相同的。

• 对于超立方体，大部分的体积位于其角落，而非中心。对于超球体，大部分的体积位于其表面，而非中心。

• 在高维空间中，多元高斯分布的大部分质量并不集中在均值（或众数）附近，而是位于一个壳层中，这个壳层从均值向尾部扩展，随着维度的增加而远离均值。这个壳层被称为典型集。

有关说明这些概念的代码示例，请访问本书的仓库：github.com/aloctavodia/BAP3。

对于我们当前的讨论，所有这些事实意味着，如果我们不明智地选择评估后验分布的位置，我们将大部分时间都花在计算对后验贡献几乎为零的值上，从而浪费宝贵的资源。因此，网格方法并不是一个明智的选择来评估后验分布，因此作为高维问题的一般方法并不十分有用。

10.3 二次方法

二次近似，也称为 Laplace 方法或正态近似，包含将后验分布近似为高斯分布。为此，我们首先找到后验分布的模型；在数值上，我们可以使用优化方法来实现。然后我们计算 Hessian 矩阵，从中可以估算标准差。如果你在想，Hessian 矩阵是一个二阶偏导数的方阵。就我们关心的内容来说，我们可以用它来获得标准差，通常是协方差矩阵的标准差。

Bambi 可以使用二次方法为我们解决贝叶斯模型。在以下代码块中，我们首先定义一个硬币投掷问题的模型，这是我们之前为网格方法定义的相同模型，然后使用 Bambi 中称为 laplace 的二次方法进行拟合：

代码 10.2

data = pd.DataFrame(data, columns=["w"]) 
priors = {"Intercept": bmb.Prior("Uniform", lower=0, upper=1)} 
model = bmb.Model("w ∼ 1", data=data, family="bernoulli", priors=priors, 
                  link="identity") 
results = model.fit(draws=4000, inference_method="laplace")

图 10.3 显示了计算得到的后验分布和精确的后验分布。请注意，Bambi 在使用此方法时还返回样本。它首先将后验分布近似为高斯分布（或多元高斯分布），然后从中采样。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file250.png

图 10.3：后验的二次近似

二次/Laplace 方法在 Bambi 中主要是出于教学目的。尽管如此，一个不错的特点是 Bambi 会考虑边界。例如，在硬币投掷问题中，我们知道解必须位于区间 [0, 1] 内。即使在背后使用高斯分布时，Bambi 也能确保这一点。Bambi 通过在一个无界的参数空间中拟合高斯分布，然后将其转换到适当的有界空间，从而实现这一点。

二次/Laplace 方法虽然本身非常有限，但可以作为更高级方法的构建块。例如，集成嵌套 拉普拉斯近似（INLA）可以非常高效地拟合多种模型。

10.4 马尔可夫方法

有一类相关的方法，统称为马尔可夫链 蒙特卡洛（MCMC）方法。这些是随机方法，只要我们能逐点计算似然函数和先验分布，就可以从真实的后验分布中获取样本。你可能记得，这正是我们在网格方法中需要的条件，但与网格方法不同，MCMC 方法能够高效地从高维空间中的高概率区域采样。

MCMC 方法根据各个区域的相对概率访问参数空间的每个区域。如果区域 A 的概率是区域 B 的两倍，那么我们从 A 中获得的样本将是从 B 中获得样本的两倍。因此，即使我们无法通过解析方式计算整个后验分布，我们仍然可以使用 MCMC 方法从中抽取样本。理论上，MCMC 将给我们来自正确分布的样本——关键在于，这一理论保证只在渐进意义上成立，也就是说，在无限数量的样本下！实际上，我们总是有有限数量的样本，因此我们需要检查这些样本是否可信。我们将要学习这一点，但不要急于求成；首先，让我们对 MCMC 方法的工作原理有些直观的了解。这将有助于我们稍后理解诊断方法。为了理解 MCMC 方法，我们将把这个方法拆分为“两个 MC 部分”：蒙特卡洛部分和马尔科夫链部分。

10.4.1 蒙特卡洛

随机数的使用解释了“蒙特卡洛”这个名字的一部分。蒙特卡洛方法是一类非常广泛的算法，通过随机抽样来计算或模拟给定的过程。蒙特卡洛是位于摩纳哥公国的一家非常著名的赌场。蒙特卡洛方法的开发者之一，Stanislaw Ulam，有一个叔叔曾在那儿赌博。Stan 的关键思想是，尽管许多问题难以精确求解或甚至无法精确表达，但它们可以通过从中抽取样本有效地进行研究。事实上，故事是这样传的：最初的动机是为了回答关于在纸牌接龙游戏中获得特定手牌的概率问题。

解决这个问题的一种方法是遵循分析组合问题。另一种方法，Stanislaw 认为，是玩几局纸牌接龙，并计算我们玩过的手牌中有多少和我们感兴趣的特定手牌相匹配！也许这对你来说听起来很显而易见，或者至少相当合理；你甚至可能用过重采样方法来解决统计问题。但请记住，这个心理实验是在大约 70 年前进行的，那时第一台实用计算机才开始被开发出来！

蒙特卡洛方法的首次应用是解决一个核物理问题，这是当时工具无法轻松应对的难题。如今，甚至个人计算机也足够强大，能够使用蒙特卡洛方法解决许多有趣的问题；因此，这些方法被广泛应用于科学、工程、工业和艺术等多个领域。使用蒙特卡洛方法计算感兴趣量的经典教学示例是对π值的数值估计。实际上，对于这个特定的计算，存在更好的方法，但它在教学中的价值依然存在。

我们可以通过以下过程估算π的值：

1. 随机地将N个点投掷到边长为 2R的正方形中。

2. 在方形内部画一个半径为R的圆，并计算圆内的点数M。

3. 计算 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/hat_Pi.png 作为比例 4https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file251.jpg。

这里有几点说明：

• 圆的面积与圆内点数(M)成正比，方形的面积与总点数(N)成正比。

• 我们知道，如果满足以下关系式，则一个点在圆内：https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file252.jpg。

• 方形的面积是(2R)²，圆的面积是πR²。因此，我们知道方形面积与圆面积的比值是π。

使用几行 Python 代码，我们可以运行这个简单的蒙特卡罗模拟并计算π，同时也能计算我们估算值相对于真实π值的相对误差。一次运行的结果如图 10.4所示。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file253.png

图 10.4：π的蒙特卡罗近似

10.4.2 马尔可夫链

马尔可夫链是一个数学对象，由一系列状态和描述如何在状态之间转换的转移概率组成。你可以自己创建一个马尔可夫链；比如，抛硬币，如果得到正面，就向右走一步，否则向左走一步。这是一个简单的一维马尔可夫链。如果一个链是马尔可夫链，那么从任何状态转移到其他状态的概率只依赖于当前状态。

作为从业者，你只需要知道马尔可夫链提供了一个框架来研究 MCMC 采样器的性质（以及其他一些有用的应用）。它们并不难理解，至少是最基本的性质并不难理解。但对于你作为模型构建者来说，深入细节并没有太大意义，因此我们不会进一步讨论。如果你有兴趣，可以查看 Blitzstein [2019]，那是一个很好的入门介绍。

最流行的 MCMC 方法可能是美特罗波利斯-哈斯廷斯算法，接下来我们将讨论它。

10.4.3 美特罗波利斯-哈斯廷斯算法

对于某些分布，例如高斯分布，我们有非常高效的算法可以从中获得样本，但对于其他分布则不一定是这样。美特罗波利斯-哈斯廷斯算法使我们能够从任何概率分布中获得样本，只要我们能够计算出与之成比例的值，因此可以忽略归一化因子。这非常有用，因为很多时候更困难的部分恰恰是计算归一化因子。这在贝叶斯统计中尤为明显，在贝叶斯统计中，计算边际似然可能是一个决定性的难点。

为了概念性地理解这种方法，我们将使用以下类比。假设我们有兴趣找到湖泊中的水量，并且找出湖泊中最深的部分。湖水非常浑浊，因此我们无法通过看水底来估算深度，湖泊又非常大，所以使用网格近似方法似乎不是一个好主意。为了制定抽样策略，我们寻求了两位好朋友的帮助：Markovia 和 Monty。经过一次富有成效的讨论，他们提出了以下算法，它需要一艘船——不用太奢华，我们甚至可以使用一只木筏和一根非常长的木棍。这比声呐便宜，而且我们已经把所有的钱都花在船上了！看看这些步骤：

1. 在湖中选择一个随机位置，并将船移到那里。

2. 使用木棍测量湖泊的深度。

3. 将船移动到另一个点并进行新的测量。

4. 通过以下方式比较这两种方法：

   1. 如果新位置比第一个位置更深，记下新位置的深度并从第 3 步重复。

   2. 如果这个位置比第一个位置更浅，我们有两个选择：接受或拒绝。接受意味着我们记录下新位置的深度并从第 3 步重复。拒绝意味着我们返回第一个位置，并再次（是的，重新！）记录第一个位置的深度值。

决定是否接受或拒绝新位置的规则称为 Metropolis-Hastings 标准，基本上说的是，我们必须按与新旧位置深度比值成正比的概率接受新位置。

如果我们按照这个迭代过程进行，不仅能得到湖泊的总体水量和最深的点，还能得到湖底整个曲率的近似值。正如你可能已经猜到的那样，在这个类比中，湖底的曲率就是后验分布，最深的点是众数。根据我们的朋友 Markovia 说，迭代次数越多，近似值就越好。事实上，理论保证在某些一般条件下，如果我们获取无限多个样本，我们最终会得到精确的答案。幸运的是，在实践中，对于许多问题，我们可以通过有限且相对较少的样本获得非常准确的近似值。

上述解释足以帮助我们概念性地理解 Metropolis-Hastings 算法。接下来的几页将提供更详细和正式的解释，以防你想深入了解。

Metropolis-Hastings 算法有以下步骤：

1. 为参数 x[i] 选择一个初始值。这可以通过随机选择或根据经验猜测来完成。

2. 通过从 Q(x[i+1]|x[i]) 中抽样，选择一个新的参数值 x[i+1]。我们可以把这一步看作是对状态 x[i] 的某种扰动。

3. 使用 Metropolis-Hastings 标准计算接受新参数值的概率：

   https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file254.jpg

4. 如果在步骤 3 中计算的概率大于从 [0, 1] 区间的均匀分布中取出的值，我们就接受新的状态；否则，我们保持在旧的状态中。

5. 我们从步骤 2 开始迭代，直到我们有足够的样本。

这里有几个需要注意的事项：

• Q 被称为提议分布。它可以是我们想要的任何分布，但我们通常选择一个容易采样的分布，比如高斯分布或均匀分布。

• 注意，Q 不是先验、似然或模型的任何部分。它是 MCMC 方法的一个组成部分，而不是模型的一部分。

• 如果 Q 是对称的，则 q(x[i]|x[i+1]) 和 q(x[i+1]|x[i]) 会相互抵消。因此，我们只需要评估比值 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file255.jpg。

• 步骤 3 和步骤 4 表明我们总是会接受转向一个更可能的状态。较不可能的参数值会按概率接受，给定新参数值 x[i+1] 和旧参数值 x[i] 之间的比率。接受提议步骤的这一标准，相较于网格方法，给我们提供了一种更高效的采样方法，同时保证了正确的采样。

• 目标分布（贝叶斯统计中的后验分布）是通过一组采样的参数值来近似的。如果我们接受，就将 x[i+1] 添加到新采样值的列表中。如果我们拒绝，就将 x[i] 添加到列表中，即使这个值是重复的。

在这个过程的最后，我们将得到一组数值。如果一切都按正确的方式完成，这些样本将是后验分布的近似。我们追踪中出现频率最高的数值将是根据后验分布最可能的值。此过程的一个优点是，分析后验分布就像操作一组数值数组一样简单，正如你在之前的章节中所做的那样。

以下代码演示了 metropolis 算法的一个非常基础的实现。它并不是为了求解一个实际问题，而是为了展示如果我们知道如何逐点计算概率密度，我们可以从概率分布中采样。请注意，以下实现没有涉及任何贝叶斯的内容；它没有先验，甚至没有数据！请记住，MCMC 方法是非常通用的算法，可以应用于广泛的问题。

metropolis 函数的第一个参数是一个 PreliZ 分布；我们假设我们不知道如何直接从这个分布中获取样本：

代码 10.3

def metropolis(dist, draws=1000): 
    """A very simple Metropolis implementation""" 
    trace = np.zeros(draws) 
    old_x = 0.5 
    old_prob = dist.pdf(old_x) 

    delta = np.random.normal(0, 0.5, draws) 
    for i in range(draws): 
        new_x = old_x + delta[i] 
        new_prob = dist.pdf(new_x) 
        acceptance = new_prob / old_prob 
        if acceptance >= np.random.random(): 
            trace[i] = new_x 
            old_x = new_x 
            old_prob = new_prob 
        else: 
            trace[i] = old_x 
    return trace

我们的简单 metropolis 算法的结果如 图 10.5 所示。黑线表示真实分布，而条形表示我们计算的样本。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file256.png

图 10.5：来自简单的 metropolis 算法的样本

算法的效率在很大程度上依赖于提议分布；如果提议的状态与当前状态相差很远，拒绝的概率就会非常高；如果提议的状态非常接近当前状态，我们就会非常缓慢地探索参数空间。在这两种情况下，我们需要的样本数量都远远多于较不极端的情况。通常，提议分布是一个多变量高斯分布，其协方差矩阵在调节阶段确定。PyMC 通过遵循一个经验法则来自适应地调节协方差，即一维高斯分布的理想接受率大约是 50%，而 n 维高斯目标分布的理想接受率大约是 23%。

MCMC 方法通常需要一些时间才能开始从目标分布中获取样本。所以，在实践中，人们会进行烧入步骤，即消除前一部分样本。进行烧入是一个实际的技巧，而不是马尔可夫理论的一部分；实际上，对于无限样本来说，这并不是必需的。因此，去除前一部分样本只是一个临时的技巧，用于在只能计算有限样本的情况下获得更好的结果。拥有理论保证或指导总是比没有它们要好，但对于任何实际问题，理解理论与实践之间的区别非常重要。记住，我们不应当因为将数学对象与这些对象的近似混淆而感到困惑。球体、高斯分布、马尔可夫链以及所有数学对象只存在于理念的柏拉图世界中，而不在我们不完美的现实世界中。

到这一点，我希望你已经对 Metropolis-Hastings 方法有了一个清晰的概念理解。你可能需要回头再读这一部分几遍；这完全没问题。主要思想简单但也很微妙。

10.4.4 哈密顿蒙特卡洛

MCMC 方法，包括 Metropolis-Hastings，提供了理论上的保证：如果我们获取足够多的样本，就能准确地近似正确的分布。然而，在实践中，可能需要比我们拥有的时间更多的时间来获取足够的样本。因此，已经提出了 Metropolis-Hastings 算法的替代方案。

许多替代方法，比如 Metropolis-Hastings 算法本身，最初是为了解决统计力学中的问题而开发的。统计力学是物理学的一个分支，研究原子和分子系统的性质，因此可以通过物理系统的类比以非常自然的方式进行解释。一个这样的修改被称为哈密顿蒙特卡洛，或者混合蒙特卡洛（HMC）。简单来说，哈密顿量是描述物理系统总能量的量。术语混合也被使用，因为它最初被构想为 Metropolis-Hastings 与分子力学的混合，分子力学是用于分子系统的广泛应用的模拟技术。

从概念上讲，我们可以将 HMC 方法视为一种 Metropolis-Hastings 方法，只不过其提议分布不是随机的。为了在不涉及数学细节的情况下对 HMC 有一个概念性的理解，我们可以再次使用湖泊和船只的类比。我们不再随机地移动船只，而是沿着湖底的曲率来移动船只。为了决定船只的移动方向，我们让一个球从当前位置开始滚到湖底。我们的球是一个非常特殊的球：它不仅是完美的球形，而且没有摩擦力，因此不会被水或泥土减速。我们扔下这个球，让它滚动一小段时间，直到我们突然停止它。然后，我们使用 Metropolis 准则接受或拒绝这个提议的步骤，就像我们在传统的 Metropolis-Hastings 方法中所做的那样。接着，整个过程会被重复多次。幸运的是，这种修改后的程序会增加接受新位置的机会，即使这些新位置相对于之前的位置远得多。

根据参数空间的曲率进行移动被证明是一种更聪明的移动方式，因为它避免了 Metropolis-Hastings 方法的主要缺点之一：高效探索样本空间需要拒绝大多数提出的步骤。相反，使用 HMC 方法，即使是参数空间中较远的点，也能获得较高的接受率，从而实现非常高效的采样方法。

让我们走出这个假想实验，回到现实世界。我们必须为这种非常巧妙的基于哈密顿量的提议付出代价。我们需要计算我们函数的梯度。梯度是导数概念在多维空间中的推广；计算函数在某一点的导数告诉我们函数在哪个方向上增大，在哪个方向上减小。我们可以使用梯度信息来模拟球体在弯曲空间中的运动；事实上，我们使用与经典物理学中模拟经典力学系统（如滚动的球体、行星系统中的轨道、分子震动）相同的运动定律和数学工具。

计算梯度使我们面临一个权衡：每一步 HMC 的计算成本比 Metropolis-Hastings 步骤高，但 HMC 的接受概率远高于 Metropolis。为了平衡这个权衡，使其有利于 HMC，我们需要调整 HMC 模型的一些参数（类似于调整 Metropolis-Hastings 采样器的提议分布宽度）。当这种调整是手动进行时，需要通过反复试验，并且还需要经验丰富的用户，这使得这个过程比我们希望的更加依赖于用户的经验，缺乏普遍适用性。幸运的是，现代的概率编程语言配备了高效的自适应哈密顿蒙特卡洛方法，比如 PyMC 中的 NUTS 采样器。这个方法已经证明在求解贝叶斯模型时非常有用和高效，且无需人工干预（或者至少可以最小化人工干预）。

哈密顿蒙特卡洛方法的一个局限是，它们仅适用于连续分布；原因是我们无法对离散分布计算梯度。PyMC 通过将 NUTS 分配给连续参数，将其他采样器分配给其他参数（例如，将 PGBART 分配给 BART 随机变量，将 Metropolis 分配给离散变量）来解决这个问题。

基于 JAX 的采样

JAX 是一个旨在提供高性能数值计算和自动微分的库，用于复杂的数学运算。PyMC 使用的是 NUTS 的 Python 版本。但你也可以使用基于 JAX 的采样器实现。根据你的模型，这些采样器可能比 PyMC 默认的 NUTS 采样器快得多。为了使用它们，我们需要为pm.sample()指定nuts_sampler参数。当前支持的选项有"nutpie"、"blackjax"和"numpyro"。这三种采样器默认不与 PyMC 一起安装，所以你需要手动安装它们。对于 CPU，nutpie 可能是最快的选择：github.com/pymc-devs/nutpie。在本书中，我们使用了 nutpie 从高斯过程（GPs）中采样——请参见第八章中的 Jupyter 笔记本。

我强烈推荐你配合这个非常酷的应用程序来阅读本节内容，作者是 Chi Feng：chi-feng.github.io/mcmc-demo。

10.5 序列蒙特卡洛

Metropolis-Hastings 和 NUTS（以及其他哈密顿蒙特卡洛变体）的一大局限是，如果后验分布具有多个峰值且这些峰值之间被非常低概率的区域分隔，这些方法可能会陷入单一模式，错过其他峰值！

许多为克服这个多重极小值问题而开发的方法基于温度化的思想。这个思想再次借鉴了统计力学的概念。一个物理系统可以占据的状态数取决于系统的温度；在 0 开尔文（最低可能温度）下，所有系统都陷入一个单一状态。而在另一极端，对于无限温度，所有可能的状态都是等概率的。通常，我们关注的是处于某个中间温度的系统。对于贝叶斯模型，通过对贝叶斯定理进行某种变化，可以直观地应用这个温度化思想。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file257.jpg

参数β被称为逆温度或温度化参数。注意，当β = 0 时，我们得到p(y|θ)^(β) = 1，因此温度化后验p(θ|y)[β]仅仅是先验p(θ)，而当β = 1 时，温度化后验就是实际的完整后验。由于从先验分布进行采样通常比从后验分布采样更容易（通过增加β的值），我们从一个较容易的分布开始采样，并逐渐将其转变为我们真正关心的更复杂的分布。

有许多方法利用了这个思想，其中之一被称为顺序蒙特卡罗（SMC）。PyMC 中实现的 SMC 方法可以总结如下（也见图 10.6）：

1. 将β初始化为 0。

2. 从温度化后验中生成N个样本S[β]。

3. 将β稍微增加一点。

4. 计算一组N的权重W。这些权重是根据新的温度化后验计算得出的。

5. 通过根据W对S[b]进行重采样，获得S[w]。

6. 运行N个 Metropolis 链，每个链都从S[w]中的不同样本开始。

7. 从步骤 3 开始重复，直到β ≥ 1。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file258.png

图 10.6：SMC 的示意图

重采样步骤通过移除概率较低的样本并用概率较高的样本替换它们来工作。Metropolis 步骤扰动这些样本，帮助探索参数空间。

温度化方法的效率很大程度上取决于β的中间值，这通常被称为降温计划。β的两个连续值之间的差异越小，两个连续的温度化后验就会越接近，因此从一个阶段到下一个阶段的过渡就会更容易。但如果步长太小，我们将需要许多中间阶段，超过某个点，这将导致浪费大量计算资源，而并未真正提高结果的准确性。

幸运的是，SMC 可以自动计算β的中间值。具体的降温计划将根据问题的难度进行调整；较难采样的分布将需要比简单分布更多的阶段。

在图 10.6 的顶部，我们有九个样本或粒子（灰色点），它们来自先验，表示为覆盖所有内容的非常宽的分布（阶段 0）。对于后续阶段，我们根据其加权后验密度重新加权前一个阶段的样本。然后，我们根据这些权重进行重采样。因此，一些粒子会被丢失并被其他样本替代，所以总数保持不变。接着，我们对样本进行突变，也就是说，我们对粒子应用一个或多个 MCMC 步骤。然后我们增加 β 并重复这一过程。当我们达到 β = 1 时，粒子（或样本）将按照后验分布进行分布。

除了β的中间值外，还有两个参数是根据前一个阶段的接受率动态计算的：每个马尔可夫链的步数和提议分布的宽度。

10.6 诊断样本

在本书中，我们使用数值方法计算几乎所有模型的后验。你在使用贝叶斯方法处理自己问题时，很可能也会遇到这种情况。由于我们使用有限数量的样本来近似后验，因此检查我们是否有有效样本是非常重要的；否则，从中得出的任何分析都会完全失真。我们可以执行几个测试，其中一些是视觉测试，其他的是定量测试。这些测试旨在发现样本的问题，但无法证明我们拥有正确的分布；它们只能提供样本似乎合理的证据。如果我们发现样本有问题，还有很多解决方案可以尝试。我们将在诊断过程中讨论这些解决方法。

为了使解释更具体，我们将使用最简模型，包含两个参数：一个全局参数 a 和一个组参数 b。仅此而已，这些模型甚至没有似然/数据！

代码 10.4

with pm.Model() as model_c: 
    a = pm.HalfNormal('a', 10) 
    b = pm.Normal('b', 0, a, shape=10) 
    idata_cm = pm.sample(tune=2000) 

with pm.Model() as model_nc: 
    a = pm.HalfNormal('a', 10) 

    b_offset = pm.Normal('b_offset', mu=0, sigma=1, shape=10) 
    b = pm.Deterministic('b', 0 + b_offset * a) 
    idata_ncm = pm.sample(tune=2000)

model_c和model_nc模型之间的区别在于，对于前者，我们直接拟合组级参数，而对于后者，我们将组级参数建模为一个平移和缩放的高斯分布。

这两个模型可能看起来对你来说过于人工化，或者只是奇怪。然而，需要注意的是，这两个模型的结构本质上与我们在第四章 4 中已经讨论过的居中和非居中参数化是相同的。

根据该章节的讨论，我们应该期待从model_nc获得比model_c更好的样本。让我们检查我们的预期是否成立。

10.7 收敛性

理论上，MCMC 方法在我们采集无限样本时是可以保证收敛的。实际上，我们需要检查所采样本是否合理且有限。通常，当我们收集到证据表明样本在某种意义上是稳定的时，我们就认为采样器已收敛。一个简单的测试方法是多次运行相同的 MCMC 模拟，并检查每次是否得到相同的结果。这也是为什么 PyMC 默认运行更多链条的原因。对于现代计算机来说，这是几乎没有成本的，因为我们有多个核心。而且，它们不会浪费资源，因为我们可以将不同链条的样本合并来计算摘要、绘图等。

有许多方法可以检查不同链条在实际应用中的等效性，无论是通过视觉方式还是正式测试。我们在这里不会深入探讨技术细节；我们只会展示一些例子，希望它们足以让你培养出解读诊断结果的直觉。

10.7.1 轨迹图

检查收敛性的一种方法是通过视觉检查链条是否相似。例如，我们可以使用 ArviZ 的 plot_trace 函数。为了更好地理解在检查这些图表时应该注意什么，我们将比较之前定义的两个模型的结果。

变量 b 是 10 维的。为了清晰和简洁，我们只展示其中一个维度。你可以随时在自己的计算机上将所有维度可视化。图 10.7 显示了许多问题。在左列中，我们有四个 KDE，每个链条一个。我们可以看到它们看起来不同。这表明每条链条采样的后验区域略有不同。在右列中，我们有轨迹本身。我们也有四条线，每条链条一条，虽然它们有些凌乱，但我们仍然可以看到，有一条链条从第一步就卡在 0 附近，直到几乎第 400 步才有所变化。我们在第 800 步左右也看到类似的情况。

当我们将图 10.7与图 10.8进行比较时，问题变得更加明显。对于后者，我们看到四条链的 KDE 看起来相互之间更为相似，而轨迹图看起来模糊得多，更像是噪声，很难看到明显的模式。我们希望曲线能够自由地漫游。当这种情况发生时，我们说我们有良好的 混合。我们这么表达是因为很难区分一条链与另一条链；它们是混合的。这是好事，因为这意味着即使我们从不同的起点启动四条（或更多）独立链条，它们也都描述了相同的分布。这并不是收敛的证明，但至少我们没有看到非收敛或糟糕混合的证据。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file259.png

图 10.7：model_c 的轨迹图

图 10.7 还有一些黑色垂直线条出现在顶部，而图 10.8 中没有这些线条。这些是发散点；本章稍后有专门的部分讨论这些问题。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file260.png

图 10.8：model_nc 的轨迹图

10.7.2 排名图

轨迹图可能很难解读，特别是当我们有多个链时，因为容易忽略一些细节。一个替代方案是使用排名图 [Vehtari et al., 2021]。为了构建某一参数的排名图，我们首先将所有链中的所有样本进行排序，并分配一个整数：这是第 0 个样本，这是第 1 个样本，这是第 2 个样本，依此类推。然后，我们将所有排名根据原始链进行分组。最后，我们绘制与链数相等的直方图。如果所有链都来自相同的分布，我们可以预期所有链都有相同数量的低排名、高排名、中等排名等。换句话说，排名的直方图应该是均匀的。

要获取排名图，我们可以调用 ArviZ 的 plot_trace 函数，并使用 kind="rank_bars" 参数。图 10.9 和图 10.10 是这类图的示例。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file261.png

图 10.9：model_c 的排名图

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file262.png

图 10.10：model_nc 的排名图

左边是我们之前展示过的相同的 KDE 图。右边是排名图。再次强调，model_nc 的结果看起来好多了；与均匀分布的偏差非常小。另一方面，我们可以从 图 10.9 中看到一些问题；例如，参数 a 的排名在 500 或更低的地方看起来非常差，参数 b 在排名大约为 2000 处也很糟糕。其他区域也存在问题。

10.7.3 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/hat_R.png (R hat)

比较独立链的一种定量方法是使用 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/hat_R.png 统计量。此测试的思路是计算链之间的方差与链内部的方差之比。理想情况下，我们应该期望值为 1。作为经验法则，值小于 1.01 也是可以接受的；较高的值则表示缺乏收敛性。我们可以使用 az.r_hat 函数计算该值（参见 表 10.1）。https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/hat_R.png 诊断量也可以通过默认的 az.summary 函数计算，或者通过 az.plot_forest（使用 r_hat=True 参数）来选择性计算。

	a	b[0]	b[1]	b[2]	b[3]	b[4]	b[5]	b[6]	b[7]	b[8]	b[9]
model_c	1.2	1.17	1.05	1.17	1.17	1.15	1.11	1.09	1.17	1.18	1.17
model_nc	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0	1.0

表 10.1：model_c 和 model_ncm 模型的 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/hat_R.png 值

大约 1.1 的值在建模初期可能是可以接受的，特别是当你只是检查一个似然是否合理，或仅仅是尝试找出你真正想要构建的模型时。此外，对于参数较多的模型来说，1.01 的阈值可能过于严格。原因是即使你确实已经达到了收敛，你仍然可能偶然得到一些超过这个阈值的 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/hat_R.png 值。例如，PyMC-BART 包括 plot_convergence 函数。该函数旨在检查 BART 随机变量的收敛性。使用 BART 模型时，每个观察值都会得到一个 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/hat_R.png，而且这个数量可能很多。因此，plot_convergence 显示了 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/hat_R.png 值的累积分布以及一个阈值，这个阈值会自动计算并考虑到观察值的数量，从而进行多重比较的校正。

图 10.11 显示了此类图的示例。在右侧，我们展示了 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/hat_R.png 的累积分布，并有一条灰色虚线表示调整后的阈值。理想情况下，整个累积曲线应位于虚线的左侧。在 图 10.11 的左侧子图中，我们展示了 有效样本量（ESS）。我们将在下一节解释 ESS。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file263.png

图 10.11：使用 pmb.plot_convergence(.) 计算的诊断图

10.8 有效样本量（ESS）

MCMC 样本可能是相关的。原因是我们使用当前的位置生成一个新位置，并在考虑旧位置的基础上接受或拒绝下一个位置。这种依赖性对于调优良好的现代方法（如哈密顿蒙特卡洛）通常较低，但有时也可能较高。我们可以使用 az.plot_autocorrelation 计算并绘制自相关图。但通常，更有用的度量是计算 有效样本量（ESS）。我们可以将这个数字看作是样本中有用的抽样次数。由于自相关，这个数字通常会低于实际的样本数。我们可以使用 az.ess 函数来计算它（见 表 10.2）。ESS 诊断也可以通过 az.summary 函数默认计算，或者通过 az.plot_forest 函数（使用 ess=True 参数）来选择性计算。

	a	b[0]	b[1]	b[2]	b[3]	b[4]	b[5]	b[6]	b[7]	b[8]	b[9]
model_cm	14	339	3893	5187	4025	5588	4448	4576	4025	4249	4973
model_ncm	2918	4100	4089	3942	3806	4171	3632	4653	3975	4092	3647

表 10.2：model_c 和 model_ncm 模型的 ESS 值

一般规则是我们至少需要一个有效样本量为 400（每条链 100 ESS）。如果得到的值低于这个数值，我们的估计不仅可能过于嘈杂，甚至像 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/hat_R.png 这样的诊断图也可能变得不可靠。

MCMC 样本的质量可能会因后验分布的不同区域而异。例如，至少对于某些问题，采样分布的中心部分可能比尾部部分更容易。因此，我们可能希望计算后验分布不同区域的 ESS 值。az.ess()返回的默认值是中心区域的 ESS，它估计分布中心的解析精度。如果你对某个参数的均值或中位数感兴趣，应该检查这个 ESS 值。如果你想报告后验区间或关注稀有事件，应检查尾部 ESS 值，它是在 5%和 95%分位数处计算的最小 ESS 值。如果你对特定分位数感兴趣，可以使用az.ess(., method='quantile')来请求这些具体的值。我们甚至可以使用az.plot_ess(., kind="quantiles")函数同时绘制多个分位数的 ESS，如图 10.12所示，针对参数a。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file264.png

图 10.12：参数a的 ESS 分位数

最后，当我们运行模型并发现 ESS 值非常低时，第一反应可能是增加样本数。有时候这就足够了。但有时即使增加 10 倍样本也不够。我们可以使用az.plot_ess(., kind="evolution")，而不是通过反复试验。这将为我们提供一个样本与 ESS 的图，如图 10.13所示。我们可以利用这些信息来估计需要多少样本才能达到给定的 ESS 值。例如，在图 10.13中，我们可以看到仅通过增加样本数，model_c中参数a的 ESS 值没有太大希望变得理想。与此相比，model_nc的 ESS 值对于大部分样本接近实际样本数。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file265.png

图 10.13：参数a的 ESS 演变

10.9 蒙特卡罗标准误差

即使我们有一个非常低的https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/hat_R.png和一个非常高的 ESS 值，MCMC 的样本仍然是有限的，因此我们在估计后验参数时会引入误差。幸运的是，我们可以估算这个误差，它被称为蒙特卡罗标准误差（MCSE）。MCSE 的估算考虑到样本之间并非真正独立。我们希望得到的结果精度受限于这个值。如果某个参数的 MCSE 为 0.2，那么报告该参数为 2.54 是没有意义的。相反，如果我们重复模拟（使用不同的随机种子），我们应该预期 68%的结果会落在 2.54±0.2 的范围内。同样，对于 95%的结果，它们应该落在 2.54±0.4 的范围内。在这里，我假设 MCSE 服从正态分布，并利用高斯分布的性质：大约 68%的值位于一个标准差以内，约 95%的值位于两个标准差以内。

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/hat_R.png、ESS 和 MCSE 是相关的。实际上，我们应该使用 ESS 作为无量纲诊断，以确保我们拥有足够的有效样本。它是无量纲的，因为无论一个参数从 0 到 1，另一个从 0 到 100，都不影响 ESS 的比较。我们可以比较它们的 ESS 值。对于 ESS，值越大越好，最低值至少为 400。如果达到最低值，我们检查是否有足够低的 https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/hat_R.png。我们还可以通过可视化的排名图或轨迹图来检查（我们还应检查发散，因为稍后我们将进行解释）。如果一切正常，我们再检查 MCSE 是否足够低，适用于我们希望报告的参数和精度。希望对于大多数问题，MCSE 都会远低于我们所需的精度。

过多的数字可能会有害

在报告文本、表格或图表时，重要的是要意识到过多的数字会让数字难以读取和理解。数字如 0.9 比 0.909297 更容易阅读，也更容易记住。此外，当报告的数字位数超过实际需要时，技术受众可能会认为你暗示着比实际存在的更高的显著性。因此，你会误导这些受众，让他们试图从本无意义的差异中寻找含义。最后，过多的数字会让你的图表、表格和图形看起来凌乱且视觉上令人不堪重负。所以，请时刻记住要考虑你的受众的数据兴趣和上下文。

10.10 发散

我们现在将探讨发散（divergences），这是一种 NUTS 特有的诊断方法，因为它基于方法的内部机制，而不是生成样本的属性。发散是一种强大且敏感的诊断方法，它表明采样器很可能已经发现了后验分布中的一个高曲率区域，该区域无法被正确探索。发散的一个优点是，它们通常出现在问题参数空间区域附近，因此我们可以利用它们来识别问题可能所在的位置。

让我们通过视觉辅助工具讨论发散：

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file266.png

图 10.14：model_c 和 model_nc 模型中选定参数的配对图

如你所见，图 10.14 显示了以下三个子图：

• 左侧子图：我们展示的是model_c模型的两个参数的散点图；即，b参数的一个维度（我们随便挑了一个 – 你可以选择不同的维度），以及a参数的对数值。我们取对数是因为a被限制为正值（它是一个尺度参数）。在采样之前，PyMC 会将所有有界的参数转换为无界的参数。对于像a这样的参数，转换方法是取对数。我们在这里也做同样的处理，因为我们想理解采样器究竟在“看到”什么。好了，我们有一个散点图，其中灰色点表示样本。看看这个参数的形状，这种形状被称为 Neal 的漏斗形状，在层次模型中很常见。黑色点表示发散，它们散布在各处，但我们可以看到，许多黑点集中在漏斗的尖端附近。这个几何形状对大多数 MCMC 方法来说是有问题的，因为很难调整采样器，使其既能在漏斗尖端和顶部获取良好的样本，又能适应这种几何形状。一种是更“球形”的区域，采样器可以上下左右自由移动，另一种是“更窄”的区域，采样器必须主要沿上下方向移动，左右方向的移动则非常有限。

• 中间的子图：我们基本上与之前的情况相同，但对于model_nc模型，现在漏斗形状更加明显。但我们没有出现发散。而且我们从前面的章节已经知道，这个模型的样本实际上更好。到底发生了什么？理解这一点的关键在于模型的定义。你会注意到，对于这个模型，b并没有被采样：b是一个确定性变量，是b_offset和a的组合，这两个参数在最后一个子图中被绘制出来。

• 右侧子图：我们有b_offset与a的关系，可以看到几何形状更“球形”。正是这一点，而不是中间的子图，才是采样器“看到”的内容。因为这种几何形状更容易采样，所以我们没有出现发散，整体的诊断结果也更好。

改变模型的参数化是去除发散的一种方法，但除非你已经知道模型的另一种参数化方法，否则找到合适的参数化可能非常耗时。一个常见且容易尝试的替代方案是改变target_accept的值，这是pm.sample的一个参数。有时候，你可能需要同时更换参数化方法和target_accept的值。那么，什么是target_accept呢？它是控制 PyMC 中 NUTS 采样器调优的一个参数。它控制提议样本的接受率，默认值为 0.8。也就是说，接受 80%的提议样本。NUTS 采样器通过自适应调整哈密顿动力学仿真中的步长来实现目标接受率。80%是一个不错的默认值，但对于某些模型，你可能想尝试更大的值，如 0.90、0.95、0.99，甚至是 0.999，如果你不想完全失去希望的话。

10.11 保持冷静，继续尝试

当诊断显示问题时，我们应该怎么做？我们应该尝试修复它们。有时候，PyMC 会提供一些关于如何修改的建议。注意这些建议，你可以节省大量的调试时间。在这里，我列出了几个常见的操作，你可以采取：

• 检查拼写错误或其他低级错误。即使是专家，也很常犯这种“低级”错误。如果你拼错了一个变量名，很可能模型根本无法运行。但有时候，错误更加微妙，你仍然得到一个语法上有效的模型，它能运行，但语义却是错误的。

• 增加样本数量。这可能对一些轻微的问题有所帮助，比如当你接近目标 ESS（或 MCSE）时，或者当^R略高于 1.01 但差距不大时。

• 从轨迹的开始部分移除一些样本。当检查轨迹图时，你可能会观察到，前几步的部分样本相比于其余部分具有明显较高或较低的值，而其他部分看起来正常。如果是这种情况，简单地移除这些前几步的样本可能就足够了。这被称为预热，这在过去是非常常见的做法。现代的采样器已经减少了对它的需求。此外，PyMC 已经会丢弃调优阶段的样本，所以这个建议现在不再像以前那样有用。

• 修改采样器参数，例如增加调优阶段的长度，或者增加 NUTS 采样器的target_accept参数。

• 转换数据。例如，对于线性回归模型，居中协变量（减去其均值）通常能加速采样器，并且减少采样问题。

• 花些时间思考你的先验分布。你不应该为了加速采样器或者去除不良诊断而调整先验。你应该用先验来编码已有的知识。但通常情况下，当你这么做时，也会让采样器的工作变得更加容易。使用像 PreliZ 和先验预测检查这样的工具，帮助你编码更好的先验分布。

• 重新参数化模型，即用一种不同但等效的方式表达模型。这并不总是容易的，但对于一些常见模型，如层级模型，你已经知道一些替代的参数化方法。

10.12 总结

在本章中，我们对一些最常用的计算后验分布的方法进行了概念性的讲解。我们特别强调了 MCMC 方法，这些方法旨在适用于任何给定的模型（或者至少广泛的模型范围），因此有时被称为通用推断引擎。这些方法是任何概率编程语言的核心，因为它们允许自动推断，从而让用户专注于迭代的模型设计和结果的解释。

我们还讨论了用于诊断样本的数值和视觉测试。没有良好的后验分布近似，贝叶斯框架的所有优势和灵活性都会消失。因此，在进行任何其他类型的分析之前，评估样本的质量是一个至关重要的步骤。

10.13 练习

1. 使用网格方法与其他先验；例如，尝试 prior = (grid <= 0.5).astype(int) 或 prior = abs(grid - 0.5)，或者尝试定义你自己的奇特先验。还可以尝试使用其他数据，例如增加数据的总量，或使数据在正反面次数上更加均匀或不均匀。

2. 在估计 π 的代码中，保持 N 固定并重新运行几次。注意到结果有所不同，因为我们使用的是随机数，但也要检查误差大致相同的顺序。尝试改变 N 点的数量，并重新运行代码。你能估算出 N 点的数量和误差之间的关系吗？为了更好的估算，你可能想修改代码，将误差作为 N 的函数进行计算。你还可以用相同的 N 运行几次代码，计算平均误差和误差的标准差。你可以使用 Matplotlib 的 plt.errorbar() 函数绘制这些结果。尝试使用一组 N 值，如 100、1000 和 10000；也就是一个数量级的差异。

3. 修改你传递给 metropolis 函数的 dist 参数；尝试使用 第一章 中的先验值。将此代码与网格方法进行比较；哪一部分需要修改，才能使用它来解决贝叶斯推断问题？

4. 将你在之前练习中的答案与 Thomas Wiecki 的代码进行比较：twiecki.github.io/blog/2015/11/10/mcmc-sampling/

5. 重温至少几个前面章节中的模型，并运行我们在本章中看到的所有诊断工具。

6. 重温所有前几章的代码，找出那些存在偏差的地方，并尽量减少这些偏差。

加入我们的社区 Discord 空间

加入我们的 Discord 社区，与志同道合的人一起学习，超过 5000 名成员在此一起成长： packt.link/bayesian

https://github.com/OpenDocCN/freelearn-ml-pt2-zh/raw/master/docs/bys-anal-py-3e/img/file1.png