PyMC 简介概述#

注意：这段文字部分基于 John Salvatier、Thomas V. Wiecki 和 Christopher Fonnesbeck 在 PeerJ CS 发表的关于 PyMC 的文章。

摘要#

概率编程允许对用户定义的概率模型进行自动贝叶斯推断。基于梯度的马尔科夫链蒙特卡罗（MCMC）采样算法，称为哈密尔顿蒙特卡罗（HMC），允许对越来越复杂的模型进行推断，但需要的梯度信息通常并不容易计算。PyMC是一个用Python编写的开源概率编程框架，利用PyTensor通过自动微分计算梯度，并即时将概率程序编译为一组计算后端之一，以提高速度。PyMC允许用Python代码进行模型规范，而不是使用领域特定语言，使其易于学习、定制和调试。本文是对这一软件包的教程式介绍，面向那些对贝叶斯统计已有一定了解的读者。

介绍#

概率编程（PP）允许以代码灵活地规范贝叶斯统计模型。PyMC是一个PP框架，具有直观且可读的强大语法，接近统计学家描述模型所使用的自然语法。它的特点是下一代马尔科夫链蒙特卡罗（MCMC）采样算法，如无转弯采样器（NUTS; Hoffman, 2014），这是一种自调节的哈密尔顿蒙特卡罗（HMC; Duane, 1987）变体。这类采样器在高维和复杂的后验分布上表现良好，允许不需要关于拟合算法的专业知识来适配许多复杂模型。HMC和NUTS利用来自似然的梯度信息，以实现比传统采样方法更快的收敛速度，特别是在较大的模型中。NUTS还有几种自调节策略，用于自适应设置哈密尔顿蒙特卡罗的可调参数，这意味着通常不需要对算法如何工作有专业知识。

Python中的概率编程带来了一系列优势，包括多平台兼容性、表现力强但简洁可读的语法、与其他科学库的易集成性，以及通过C、C++、Fortran或Cython的可扩展性。这些特性使得编写和使用自定义统计分布、采样器和变换函数相对简单，这些在贝叶斯分析中是必要的。

虽然PyMC的大部分用户界面特性都是用纯Python编写的，但它利用PyTensor（Theano项目的一个分支）透明地将模型转码为C并编译为机器代码，从而提升性能。PyTensor是一个允许使用称为张量的广义向量数据结构定义表达式的库，这些张量与流行的NumPy ndarray数据结构紧密集成，并同样允许广播和高级索引，就像NumPy数组一样。PyTensor还自动优化似然的计算图以提高速度，并允许编译为一组计算后端，包括Jax和Numba。

在这里，我们介绍如何使用PyMC解决一般贝叶斯统计推断和预测问题的基础知识。我们将首先展示如何使用PyMC的基础知识，通过一个简单的例子进行激励：安装、数据创建、模型定义、模型拟合和后验分析。接下来，我们将涵盖两个案例研究，并利用它们展示如何定义和拟合更复杂的模型。最后，我们将讨论其他几个有用的特性：自定义分布和任意确定性节点。

安装#

运行PyMC需要相对较新的Python解释器，最好是版本3.8或更高。对于macOS、Linux和Windows，最简单的获取完整Python安装的方法是下载并安装免费的Anaconda Python Distribution（由ContinuumIO提供）或开源的Miniforge。

安装Python后，请按照PyMC文档网站上的安装指南进行操作。

PyMC在宽松的Apache许可证2.0下发布。在GitHub网站上，用户也可以报告错误和其他问题，以及为项目贡献文档或代码，我们对此表示积极鼓励。

一个激励性的例子：线性回归#

为了引入模型定义、拟合和后验分析，我们首先考虑一个简单的贝叶斯线性回归模型，该模型对于参数具有正态分布的先验。我们感兴趣的是预测结果 $Y$ ，它被视为具有期望值 $μ$ 的正态分布观测值， $μ$ 是两个预测变量 $X_{1}$ 和 $X_{2}$ 的线性函数：

\begin{array}{r} \begin{aligned} Y & \sim N (μ, σ^{2}) \\ μ & = α + β_{1} X_{1} + β_{2} X_{2} \end{aligned} \end{array}

其中 $α$ 是截距， $β_{i}$ 是协变量 $X_{i}$ 的系数，而 $σ$ 表示观测误差。由于我们正在构建一个贝叶斯模型，我们必须为模型中的未知变量分配一个先验分布。我们选择零均值的正态先验，方差为 100，适用于两个回归系数，这对应于关于真实参数值的弱信息。我们选择半正态分布（在零处限制的正态分布）作为 $σ$ 的先验。

\begin{array}{r} \begin{aligned} α & \sim N (0, 100) \\ β_{i} & \sim N (0, 100) \\ σ & \sim | N (0, 1) | \end{aligned} \end{array}

生成数据#

我们可以使用仅仅 NumPy 的 random 模块从这个模型中模拟一些人工数据，然后使用 PyMC 尝试恢复相应的参数。我们故意生成数据，以便与 PyMC 模型结构紧密对应。

import arviz as az
import pandas as pd
import matplotlib.pyplot as plt
import numpy as np

%config InlineBackend.figure_format = 'retina'
# 初始化随机数生成器
RANDOM_SEED = 8927
rng = np.random.default_rng(RANDOM_SEED)
az.style.use("arviz-darkgrid")

# 真实参数值
alpha, sigma = 1, 1
beta = [1, 2.5]

# 数据集大小
size = 100

# 预测变量
X1 = np.random.randn(size)
X2 = np.random.randn(size) * 0.2

# 模拟结果变量
Y = alpha + beta[0] * X1 + beta[1] * X2 + rng.normal(size=size) * sigma

这里是模拟数据的样子。我们使用绘图库 matplotlib 中的 pylab 模块。

fig, axes = plt.subplots(1, 2, sharex=True, figsize=(10, 4))
axes[0].scatter(X1, Y, alpha=0.6)
axes[1].scatter(X2, Y, alpha=0.6)
axes[0].set_ylabel("Y")
axes[0].set_xlabel("X1")
axes[1].set_xlabel("X2");

../../_images/e65f9c582f1bb96a1507903369d4717997f318cbdb802321b2211fda7880108f.png

后验分析#

PyMC 的绘图和诊断功能现在由一个专门的、与平台无关的包 Arviz 处理。可以使用 plot_trace 创建简单的后验图。

az.plot_trace(idata, combined=True);

../../_images/cf29ae046fa23c4ea2b5d310b3cccd94f9850addcedcd2ad9ff684a03f0b944d.png

左列由每个随机变量的边际后验的平滑直方图（使用核密度估计）组成，而右列包含按顺序绘制的马尔可夫链样本。beta变量是向量值的，因此产生两个密度图和两个轨迹图，对应于两个预测系数。

此外，summary函数提供了常见后验统计量的文本输出：

az.summary(idata, round_to=2)

	mean	sd	hdi_3%	hdi_97%	mcse_mean	mcse_sd	ess_bulk	ess_tail	r_hat
alpha	1.16	0.10	0.97	1.35	0.00	0.00	6521.35	3554.23	1.0
beta[0]	0.94	0.12	0.71	1.15	0.00	0.00	6085.71	3356.17	1.0
beta[1]	2.98	0.55	1.94	3.99	0.01	0.01	5458.07	3139.07	1.0
sigma	1.01	0.07	0.87	1.13	0.00	0.00	6554.05	3380.64	1.0

案例研究 1：听障儿童的教育结果#

作为一个激励的例子，我们将使用一个听障儿童教育结果的数据集。在这里，我们感兴趣的是确定与更好或更差学习结果相关的因素。

数据#

该匿名数据集取自听力与口语语言数据 Repository (LSL-DR)，这是一个国际数据存储库，追踪听力损失儿童的人口统计信息和长期结果，这些儿童参与了支持听力和口语语言发展的项目。研究人员希望发现与这些项目中教育成果改善相关的因素。

可用的预测因素包括：

性别（male）
家庭中兄弟姐妹的数量（siblings）
家庭参与指数（family_inv）
主要家庭语言是否为英语以外的语言（non_english）
是否有先前的残疾（prev_disab）
非白人种族（non_white）
检测时的年龄（以月为单位，age_test）
听力损失是否不严重（non_severe_hl）
受试者的母亲是否获得高中文凭或以上（mother_hs）
听力障碍是否在3个月大时被识别（early_ident）。

结果变量是在几个学习领域中的标准化测试分数。

test_scores = pd.read_csv(pm.get_data("test_scores.csv"), index_col=0)
test_scores.head()

	Unnamed: 0	score	male	siblings	family_inv	non_english	prev_disab	age_test	non_severe_hl	mother_hs	early_ident	non_white
0	0	40	0	2.0	2.0	False	NaN	55	1.0	NaN	False	False
1	1	31	1	0.0	NaN	False	0.0	53	0.0	0.0	False	False
2	2	83	1	1.0	1.0	True	0.0	52	1.0	NaN	False	True
3	3	75	0	3.0	NaN	False	0.0	55	0.0	1.0	False	False
4	5	62	0	0.0	4.0	False	1.0	50	0.0	NaN	False	False

test_scores["score"].hist();

../../_images/0a5c4058a9a84ce8d699abbfcd22b6d95583c851b5dd016b2ee4028d22c32fe8.png

# 通常情况下，丢弃缺失值是一个非常糟糕的主意，但为了简化，我们在这里这样做。
X = test_scores.dropna().astype(float)
y = X.pop("score")

# 标准化功能
X -= X.mean()
X /= X.std()

N, D = X.shape

模型#

这是一个比第一个回归示例更现实的问题，因为我们现在处理的是一个 多元回归 模型。然而，尽管在 LSL-DR 数据集中有多个潜在的预测变量，但在构建有效的统计模型时，事前确定哪些变量是相关的仍然是困难的。有许多进行变量选择的方法，但一种流行的自动化方法是 正则化，通过正则化（一种惩罚形式）将无效的协变量缩小到接近零，如果它们对预测结果没有贡献。

你可能听说过正则化，这在机器学习或经典统计应用中是常见的方法，像 lasso 或 ridge 回归通过对回归参数的大小施加惩罚将参数缩小到零。在贝叶斯背景下，我们对回归系数应用适当的先验分布。其中一种先验是 分层正则化的马蹄形分布，它使用两种正则化策略，一种是全局的，另一组是局部参数，每个系数一个。要实现这一点，关键在于选择一种长尾分布作为收缩先验，这允许某些系数不为零，同时将其余的推向零。

每个回归系数 $β_{i}$ 的马蹄形先验看起来像这样：

β_{i} \sim N (0, τ^{2} \cdot {\tilde{λ}}_{i}^{2})

其中 $σ$ 是对误差标准差的先验，这也将用于模型的似然性。在这里， $τ$ 是全局收缩参数， ${\tilde{λ}}_{i}$ 是局部收缩参数。我们先从全局开始：对于 $τ$ 的先验，我们将使用一个 Half-StudentT 分布，这是一种合理的选择，因为它是重尾的。

τ \sim {Half-StudentT}_{2} (\frac{D_{0}}{D - D_{0}} \cdot \frac{σ}{\sqrt{N}}) .

一个问题是，先验的参数化要求预先指定一个值 $D_{0}$ ，该值表示非零系数的真实数量。幸运的是，对于这个值的合理猜测是所需的唯一条件，而且它只需要在真实数量的数量级内。我们将使用预测变量数量的一半作为我们的猜测：

D0 = int(D / 2)

与此同时，局部收缩参数由以下比率定义：

{\tilde{λ}}_{i}^{2} = \frac{c^{2} λ_{i}^{2}}{c^{2} + τ^{2} λ_{i}^{2}} .

为了完成这个模型定义，我们需要对 $λ_{i}$ 和 $c$ 设定先验；与全局收缩类似，我们对 $λ_{i}$ 使用长尾的 ${Half-StudentT}_{5} (1)$ 作为先验。我们需要 $c^{2}$ 严格为正，但不一定是长尾分布，因此 $c^{2}$ 的逆伽马分布先验 $c^{2} \sim InverseGamma (1, 1)$ 符合要求。

最后，为了使NUTS采样器能更高效地采样 $β_{i}$ ，我们将其重新参数化如下：

\begin{array}{r} \begin{aligned} z_{i} & \sim N (0, 1), \\ β_{i} & = z_{i} \cdot τ \cdot \tilde{λ_{i}} . \end{aligned} \end{array}

在实际操作中，你会经常遇到这种重新参数化。

模型规范#

在PyMC中指定模型与其统计规范相似。该模型采用了几个新的分布：HalfStudentT分布用于 $τ$ 和 $λ$ 先验，InverseGamma分布用于 $c^{2}$ 变量。

在PyMC中，具有完全正值先验的变量，如InverseGamma，会使用对数变换。这使得采样更为稳健。在后台，将一个无约束空间中的变量（命名为<variable-name>_log）添加到模型中用于采样。具有双侧约束的先验变量，如Beta或Uniform，也会被转化为无约束形式，但使用对数赔率变换。

我们还将利用PyMC和ArviZ中的命名维度，通过将输入变量名称作为称为“预测器”的坐标传递到模型中。这将允许我们将这个名称向量作为向量值参数的shape整数参数的替代。这一模型随后将把适当的名称与它所估计的每个潜在参数相关联。这虽然需要更多的设置工作，但在我们处理模型输出时将带来丰厚的回报。

让我们在PyMC中编码这个模型：

import pytensor.tensor as at

with pm.Model(coords={"predictors": X.columns.values}) as test_score_model:
    # 先前误差标准差
    sigma = pm.HalfNormal("sigma", 25)

    # 全局收缩先验
    tau = pm.HalfStudentT("tau", 2, D0 / (D - D0) * sigma / np.sqrt(N))
    # 局部收缩先验
    lam = pm.HalfStudentT("lam", 5, dims="predictors")
    c2 = pm.InverseGamma("c2", 1, 1)
    z = pm.Normal("z", 0.0, 1.0, dims="predictors")
    # 收缩系数
    beta = pm.Deterministic(
        "beta", z * tau * lam * at.sqrt(c2 / (c2 + tau**2 * lam**2)), dims="predictors"
    )
    # 截距无收缩
    beta0 = pm.Normal("beta0", 100, 25.0)

    scores = pm.Normal("scores", beta0 + at.dot(X.values, beta), sigma, observed=y.values)

注意，我们已经将 beta 的计算封装在 Deterministic PyMC 类中。您可以在下面更详细地阅读关于这个的内容，但这确保了这个确定性变量的值在样本轨迹中得以保留。

还要注意，我们在上下文管理器的第一次出现中声明了 Model 名称 test_score_model，而不是像第一个示例那样将其拆分成两行。

一旦模型完成，我们可以使用 GraphViz 查看其结构，这样可以绘制模型图。这有助于确保您编写的模型中的关系是正确的，因为很容易出现编码错误。

pm.model_to_graphviz(test_score_model)

../../_images/8de474e1f2035cbe504b1b5820d3b08e798723c76b3a94afc6e6a4bf3527e87e.svg

在继续之前，让我们先看看模型在未见任何数据之前的表现。我们可以进行先验预测采样，从模型中生成模拟数据。然后，让我们将这些模拟结果与数据集中实际的测试分数进行比较。

with test_score_model:
    prior_samples = pm.sample_prior_predictive(100)

Sampling: [beta0, c2, lam, scores, sigma, tau, z]

az.plot_dist(
    test_scores["score"].values,
    kind="hist",
    color="C1",
    hist_kwargs=dict(alpha=0.6),
    label="observed",
)
az.plot_dist(
    prior_samples.prior_predictive["scores"],
    kind="hist",
    hist_kwargs=dict(alpha=0.6),
    label="simulated",
)
plt.xticks(rotation=45);

../../_images/b3ec277bb08df2c650059692c23859250211a1810d3a1cb6912152a3faabef77.png

我们如何知道这是否合理？这需要对问题的领域知识。在这里，我们试图预测考试分数的结果。如果我们的模型预测的值在几千或许多负值之间，并且排除了合理的分数，那么我们就错误地设定了模型。您可以看到，模拟数据的分布支持完全重叠于观测到的分数分布的支持；这是一个好迹象！有一些负值和一些可能过大而不合理的值，但没有什么好担心的。

模型拟合#

现在进入简单的部分：PyMC 的“推断按钮”是调用 sample。我们将让这个模型调整的时间比默认值（1000 次迭代）稍长一些。这为 NUTS 采样器提供了更多的时间来充分调整自己。

with test_score_model:
    idata = pm.sample(1000, tune=2000, random_seed=42)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [sigma, tau, lam, c2, z, beta0]

100.00% [12000/12000 00:04<00:00 Sampling 4 chains, 0 divergences]

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

请注意，我们在这里有一些关于发散的警告。这些是NUTS无法在后验分布中进行有效移动的样本，因此生成的点可能不是后验的代表性样本。在这个例子中这样的样本不多，所以不必过于担心，但让我们按照建议将 target_accept 从默认值0.9增加到0.99。

with test_score_model:
    idata = pm.sample(1000, tune=2000, random_seed=42, target_accept=0.99)

Auto-assigning NUTS sampler...
Initializing NUTS using jitter+adapt_diag...
Multiprocess sampling (4 chains in 4 jobs)
NUTS: [sigma, tau, lam, c2, z, beta0]

100.00% [12000/12000 00:13<00:00 Sampling 4 chains, 0 divergences]

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

由于目标接受率更高，算法在进行跃迁步骤时变得更加保守，使得步骤变小。我们为此付出的代价是每个样本的完成时间变长。然而，警告现在已经消失，我们得到了一个干净的后验样本！

模型检测#

模型检测的一个简单第一步是通过查看单变量潜在参数的轨迹图来对我们的样本进行目视检查，以检查明显的问题。这些名称可以在 var_names 参数中传递给 plot_trace。

az.plot_trace(idata, var_names=["tau", "sigma", "c2"]);

../../_images/d0ecd18828a92325732d327aaa110c645c52a8e9bb32863522d2829d23564eb6.png

这些看起来还行吗？左侧每个参数的密度图看起来都与其他参数相似，这意味着它们已经收敛到相同的后验分布（无论它是否正确）。右侧的轨迹图的均匀性也是一个好兆头；采样值的时间序列没有趋势或模式。注意到 c2 和 tau 偶尔会采样到极端值，但这是重尾分布的预期结果。

下一个简单的模型检查步骤是查看NUTS采样器是否表现如预期。能量图是一种检查NUTS算法是否能够充分探索后验分布的方法。如果没有，可能会导致后验估计有偏，因为后验的某些部分没有被充分访问。该图显示了两个密度估计：一个是采样运行的边际能量分布，另一个是步骤之间能量转移的分布。这一切都有点抽象，但我们要寻找的只是这些分布相似。我们的结果看起来还不错。

az.plot_energy(idata);

../../_images/f5ab9d910092d5810003453555443b6ef4b1c8b8a3f916cdb7349389b6d0dc4d.png

最终，我们关心的是 beta 的估计值，即预测系数的集合。将 beta 传递给 plot_trace 会生成一个非常拥挤的图，因此我们将使用 plot_forest，它专门设计用于处理向量值参数。

az.plot_forest(idata, var_names=["beta"], combined=True, hdi_prob=0.95, r_hat=True);

../../_images/9d8cbebaeff4ff0792fd03b4f04913bf02391a1095daeb9998bd8e3e5ce36f51.png

系数的后验分布揭示了一些在预测测试分数时重要的因素。家庭参与（family_inv）值很大且为负，这意味着较高的分数（与较差的参与相关）导致测试分数显著下降。另一方面，早期识别听力障碍是正向的，这意味着及早发现问题可以在未来带来更好的教育结果，这也是直观的。值得注意的是，其他变量，特别是性别（male）、测试时年龄（age_test）和母亲的教育水平（mother_hs）几乎都被压缩到零。

案例研究 2：煤矿灾难#

考虑从1851年至1962年在英国记录的煤矿灾难的时间序列（Jarrett, 1979）。灾难的数量被认为受到这段时间安全法规变化的影响。不幸的是，我们还有一对缺失数据的年份，这些缺失值在pandas Series中被标识为nan。这些缺失值将由PyMC自动填补。

接下来，我们将为这个系列建立一个模型，尝试估计何时发生了变化。同时，我们将学习如何处理缺失数据，使用多个采样器以及从离散随机变量中进行采样。

# 格式化：关闭
disaster_data = pd.Series(
    [4, 5, 4, 0, 1, 4, 3, 4, 0, 6, 3, 3, 4, 0, 2, 6,
    3, 3, 5, 4, 5, 3, 1, 4, 4, 1, 5, 5, 3, 4, 2, 5,
    2, 2, 3, 4, 2, 1, 3, np.nan, 2, 1, 1, 1, 1, 3, 0, 0,
    1, 0, 1, 1, 0, 0, 3, 1, 0, 3, 2, 2, 0, 1, 1, 1,
    0, 1, 0, 1, 0, 0, 0, 2, 1, 0, 0, 0, 1, 1, 0, 2,
    3, 3, 1, np.nan, 2, 1, 1, 1, 1, 2, 4, 2, 0, 0, 1, 4,
    0, 0, 0, 1, 0, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1]
)
# fmt: 开启
years = np.arange(1851, 1962)

plt.plot(years, disaster_data, "o", markersize=8, alpha=0.4)
plt.ylabel("Disaster count")
plt.xlabel("Year");

../../_images/f548574c1a112c12b2c26444fc0e09af6f840e5e69fe33309378751672cae8fc.png

在时间序列中灾害的发生被认为遵循一个具有大速率参数的泊松过程，在时间序列的早期部分，而在后期部分则遵循一个具有较小速率的过程。我们有兴趣于确定序列中的变化点，这可能与矿业安全法规的变化有关。

在我们的模型中，

\begin{array}{r} \begin{aligned} D_{t} & \sim Pois (r_{t}), r_{t} = {\begin{cases} e, & 如果 t \leq s \\ l, & 如果 t > s \end{cases} \\ s & \sim Unif (t_{l}, t_{h}) \\ e & \sim exp (1) \\ l & \sim exp (1) \end{aligned} \end{array}

参数的定义如下：

$D_{t}$ : 第 $t$ 年的灾害数量
$r_{t}$ : 第 $t$ 年灾害的泊松分布的速率参数。
$s$ : 速率参数变化的年份（切换点）。
$e$ : 切换点 $s$ 之前的速率参数。
$l$ : 切换点 $s$ 之后的速率参数。
$t_{l}$ , $t_{h}$ : 年 $t$ 的下界和上界。

该模型的构建与我们之前的模型类似。主要的不同之处在于引入了具有泊松和离散均匀先验的离散变量，以及具有新颖形式的确定性随机变量 rate。

with pm.Model() as disaster_model:
    switchpoint = pm.DiscreteUniform("switchpoint", lower=years.min(), upper=years.max())

    # 切换前后的灾难发生率先验
    early_rate = pm.Exponential("early_rate", 1.0)
    late_rate = pm.Exponential("late_rate", 1.0)

    # 为当前年份前后的年份分配适当的泊松速率
    rate = pm.math.switch(switchpoint >= years, early_rate, late_rate)

    disasters = pm.Poisson("disasters", rate, observed=disaster_data)

/Users/jthompson1/miniconda3/lib/python3.10/site-packages/pymc/model/core.py:1365: ImputationWarning: Data in disasters contains missing values and will be automatically imputed from the sampling distribution.
  warnings.warn(impute_message, ImputationWarning)

速率随机变量的逻辑，

rate = switch(switchpoint >= year, early_rate, late_rate)

是通过使用 switch 实现的，该函数的工作方式类似于 if 语句。它使用第一个参数在接下来的两个参数之间进行切换。

缺失值通过在创建观察到的随机变量时将一个 NumPy MaskedArray 或一个包含 NaN 值的 DataFrame 传递给 observed 参数来透明地处理。在幕后，另一个随机变量 disasters.missing_values 被创建以建模缺失值。

不幸的是，由于它们是离散变量，因此没有有意义的梯度，我们无法使用NUTS来对switchpoint或缺失的灾害观察值进行采样。相反，我们将使用Metropolis步骤方法进行采样，该方法实现了自适应Metropolis-Hastings，因为它旨在处理离散值。PyMC会自动分配正确的采样算法。

with disaster_model:
    idata = pm.sample(10000)

Multiprocess sampling (4 chains in 4 jobs)
CompoundStep
>CompoundStep
>>Metropolis: [switchpoint]
>>Metropolis: [disasters_unobserved]
>NUTS: [early_rate, late_rate]

100.00% [44000/44000 00:05<00:00 Sampling 4 chains, 0 divergences]

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

在下面的轨迹图中，我们可以看到大约有10年的时间跨度是安全性显著变化的合理范围，但有5年的时间跨度包含了大部分的概率质量。分布看起来不平滑是因为年份切换点与可能性之间存在跳跃关系；这种锯齿状并不是由于采样误差造成的。

axes_arr = az.plot_trace(idata)
plt.draw()
for ax in axes_arr.flatten():
    if ax.get_title() == "switchpoint":
        labels = [label.get_text() for label in ax.get_xticklabels()]
        ax.set_xticklabels(labels, rotation=45, ha="right")
        break
plt.draw()

../../_images/9c66014cc14268f2c84bfb12205d98195a1bdee6ebdbe3d7afa83763f8045745.png

请注意，rate随机变量在跟踪中没有出现。这在这种情况下是可以的，因为它本身并不重要。请记住，在之前的例子中，我们通过将其包装在一个Deterministic类中并给它命名来跟踪变量。

下图显示了作为橙色垂直线的切换点，以及作为半透明带的最高后验密度（HPD）。虚线黑线显示了事故率。

plt.figure(figsize=(10, 8))
plt.plot(years, disaster_data, ".", alpha=0.6)
plt.ylabel("Number of accidents", fontsize=16)
plt.xlabel("Year", fontsize=16)

trace = idata.posterior.stack(draws=("chain", "draw"))

plt.vlines(trace["switchpoint"].mean(), disaster_data.min(), disaster_data.max(), color="C1")
average_disasters = np.zeros_like(disaster_data, dtype="float")
for i, year in enumerate(years):
    idx = year < trace["switchpoint"]
    average_disasters[i] = np.mean(np.where(idx, trace["early_rate"], trace["late_rate"]))

sp_hpd = az.hdi(idata, var_names=["switchpoint"])["switchpoint"].values
plt.fill_betweenx(
    y=[disaster_data.min(), disaster_data.max()],
    x1=sp_hpd[0],
    x2=sp_hpd[1],
    alpha=0.5,
    color="C1",
)
plt.plot(years, average_disasters, "k--", lw=2);

../../_images/b389cea58634d366755e194c206b9574ffc2a2aeb1ca36c4cc736b7e9d3521ec.png

任意确定性#

由于依赖于PyTensor，PyMC提供了许多数学函数和运算符，用于将随机变量转换为新的随机变量。然而，PyTensor中的函数库并不全面，因此PyTensor和PyMC提供了在纯Python中创建任意函数的功能，并将这些函数包含在PyMC模型中。这可以通过as_op函数修饰符来实现。

PyTensor需要知道函数的输入和输出类型，这些类型通过as_op中的itypes来指定输入类型，otypes来指定输出类型。

from pytensor.compile.ops import as_op


@as_op(itypes=[at.lscalar], otypes=[at.lscalar])
def crazy_modulo3(value):
    if value > 0:
        return value % 3
    else:
        return (-value + 1) % 3


with pm.Model() as model_deterministic:
    a = pm.Poisson("a", 1)
    b = crazy_modulo3(a)

这种方法一个重要的缺点是，pytensor无法检查这些函数以计算哈密顿基采样器所需的梯度。因此，对于使用此类运算符的模型，无法使用HMC或NUTS采样器。然而，如果我们继承自 Op 而不是使用 as_op，则可以添加梯度。PyMC示例集中包含了as_op使用的更复杂示例。

任意分布#

类似地，PyMC中的统计分布库并不详尽，但PyMC允许创建用户定义的函数以实现任意概率分布。对于简单的统计分布，DensityDist类接受一个参数，即计算对数概率 $l o g (p (x))$ 的任何函数。该函数可以在计算中使用其他随机变量。以下是一个示例，灵感来自Jake Vanderplas关于在进行线性回归时使用何种先验的博客文章（Vanderplas, 2014）。

import pytensor.tensor as at

with pm.Model() as model:
    alpha = pm.Uniform('intercept', -100, 100)
    
    # 创建自定义密度
    beta = pm.DensityDist('beta', logp=lambda value: -1.5 * at.log(1 + value**2))
    eps = pm.DensityDist('eps', logp=lambda value: -at.log(at.abs_(value)))
    
    # 创建似然
    like = pm.Normal('y_est', mu=alpha + beta * X, sigma=eps, observed=Y)

对于更复杂的分布，可以创建一个{class}~pymc.Continuous或Discrete的子类，并提供自定义的logp函数，按需实现。这就是PyMC中内置分布的指定方式。例如，心理学和天体物理学等领域对于特定过程具有复杂的似然函数，这可能需要数值近似。

下面展示了将上面的beta变量作为Continuous子类实现的代码，以及一个相关的{class}~pytensor.RandomVariable`对象，其中的一个实例成为分布的一个属性。

class BetaRV(at.random.op.RandomVariable):
    name = "beta"
    ndim_supp = 0
    ndims_params = []
    dtype = "floatX"

    @classmethod
    def rng_fn(cls, rng, size):
        raise NotImplementedError("Cannot sample from beta variable")


beta = BetaRV()

class Beta(pm.Continuous):
    rv_op = beta

    @classmethod
    def dist(cls, mu=0, **kwargs):
        mu = at.as_tensor_variable(mu)
        return super().dist([mu], **kwargs)

    def logp(self, value):
        mu = self.mu
        return beta_logp(value - mu)


def beta_logp(value):
    return -1.5 * at.log(1 + (value) ** 2)


with pm.Model() as model:
    beta = Beta("beta", mu=0)

如果您的 logp 不能用 PyTensor 表达，您可以按如下方式用 as_op 装饰该函数：@as_op(itypes=[at.dscalar], otypes=[at.dscalar])。请注意，这将创建一个黑箱 Python 函数，其速度会慢很多，并且不会提供 NUTS 等所需的梯度。

讨论#

概率编程是统计学习中的一种新兴范例，其中贝叶斯建模是一个重要的子学科。概率编程的标志性特征——将变量指定为概率分布并根据其他变量及观察结果对变量进行条件化——使其成为在各种环境和不同模型复杂度下构建模型的强大工具。随着概率编程的兴起，贝叶斯模型的拟合方法经历了一次创新的爆发，相较于现有的MCMC方法表现出显著的改进。然而，尽管这种扩展存在，能够跟上方法学创新的可用软件包仍然很少，更少的软件包允许非专业用户实现模型。

PyMC 为定量研究人员提供了一个概率编程平台，灵活且简洁地实现统计模型。丰富的统计分布库和若干预定义的拟合算法使用户可以专注于当前的科学问题，而不是贝叶斯建模的实现细节。选择 Python 作为开发语言，而不是特定领域的语言，意味着 PyMC 用户能够以交互方式构建模型、检查模型对象，并调试或分析他们的工作，使用一种易于学习的动态高级编程语言。PyMC 的模块化面向对象设计意味着添加新的拟合算法或其他功能非常简单。此外，PyMC 还具备许多其他软件包所没有的功能，最显著的是基于哈密顿的采样器以及仅由 Stan 提供的受限随机变量的自动变换。然而，与 Stan 不同的是，PyMC 支持离散变量以及非梯度基础采样算法，如梅特罗波利斯-哈斯廷斯（Metropolis-Hastings）和切片采样（Slice sampling）。

PyMC 的开发是一个持续的努力，未来版本中计划加入几个新功能。特别是，变分推断技术通常比 MCMC 采样更高效，但这以泛化能力为代价。然而，最近已经开发出了黑箱变分推断算法，例如自动微分变分推断（ADVI; Kucukelbir 等，2017）。该算法计划加入到 PyMC 中。作为一个开源科学计算工具包，我们鼓励研究人员为贝叶斯模型开发新的拟合算法，以便在 PyMC 中提供参考实现。由于采样器可以用纯 Python 代码编写，因此可以通用地实现，使它们能够在任意 PyMC 模型上工作，从而让作者接触到更大的受众来使用他们的方法。

参考文献#

Bastien, F., Lamblin, P., Pascanu, R., Bergstra, J., Goodfellow, I., Bergeron, A., Bouchard, N., Warde-Farley, D., 和 Bengio, Y. (2012) “Theano: 新特性和速度提升”。NIPS 2012 深度学习研讨会。

Bergstra, J., Breuleux, O., Bastien, F., Lamblin, P., Pascanu, R., Desjardins, G., Turian, J., Warde-Farley, D., 和 Bengio, Y. (2010) “Theano: 一个 CPU 和 GPU 数学表达式编译器”。2010 年 Python 科学计算会议 (SciPy) 论文集。2010 年 6 月 30 日至 7 月 3 日，德克萨斯州奥斯汀。

Duane, S., Kennedy, A. D., Pendleton, B. J., 和 Roweth, D. (1987) “混合蒙特卡洛”，物理快报，第 195 卷，第 216-222 页。

Hoffman, M. D., & Gelman, A. (2014)。无转弯采样器：自适应设置哈密顿蒙特卡洛中的路径长度。《机器学习研究期刊》，30。

Jarrett, R.G. 关于煤矿灾难间隔的说明。《生物计量学》，66:191–193，1979。

Kucukelbir A, Dustin Tran, Ranganath R, Gelman A, 和 Blei DM. 自动微分变分推断，《机器学习研究期刊》。18 ，第 430-474 页。

Lunn, D.J., Thomas, A., Best, N., 和 Spiegelhalter, D. (2000) WinBUGS – 一种贝叶斯建模框架：概念、结构和可扩展性。《统计与计算》，10:325–337。

Neal, R.M. 切片采样。《统计年鉴》。 (2003). doi:10.2307/3448413。
Patil, A., D. Huard 和 C.J. Fonnesbeck. (2010) PyMC：在 Python 中的贝叶斯随机建模。《统计软件期刊》，35(4)，第 1-81 页。

Piironen, J., & Vehtari, A. (2017). 稀疏性信息和其他收缩先验中的正则化。电子统计学期刊，11(2)，5018-5051。

Stan 开发团队. (2014). Stan：概率和采样的 C++ 库，版本 2.5.0。

Vanderplas, Jake. “频率主义与贝叶斯主义 IV：如何在 Python 中成为一名贝叶斯。 “Pythonic Perambulations. N.p., 2014 年 6 月 14 日。Web. 2015 年 5 月 27 日。

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

Last updated: Tue Sep 26 2023

Python implementation: CPython
Python version       : 3.10.9
IPython version      : 8.10.0

xarray: 2023.2.0

arviz     : 0.16.1
numpy     : 1.22.3
pymc      : 5.8.2
pytensor  : 2.16.3
pandas    : 1.4.2
matplotlib: 3.7.2

Watermark: 2.4.3