pymc和贝叶斯模型编程（2）

pymc中的变分推断

pymc和贝叶斯模型编程（2）。

简介和安装

简介

PyMC是一个 Python 概率编程库，允许用户使用简单的 Python API 构建贝叶斯模型，并使用马尔可夫链蒙特卡罗 (MCMC) 方法对其进行拟合。

PyMC 致力于使贝叶斯建模尽可能简单、轻松，让用户能够专注于他们的问题而不是方法。

具有如下特性：

• 现代：包括最先进的推理算法，包括 MCMC (NUTS) 和变分推断 (ADVI)。
• 用户友好：使用友好的 Python 语法编写模型。从许多示例笔记本中学习贝叶斯建模。
• 快速：使用PyTensor作为计算后端，通过 C、Numba 或 JAX 进行编译，在 GPU 上运行模型，并从复杂的图形优化中受益。
• 拓展性好：包括概率分布、高斯过程、ABC、SMC 等等。它与用于可视化和诊断的ArviZ以及用于高级混合效果模型的Bambi完美集成。
• 以社区为中心：在讨论中提出问题、加入MeetUp 活动、在Twitter上关注我们并开始贡献。

PyMC和PyMC3？

PyMC3已重命名为PyMC，PyMC3 基于Theano计算 ，PyMC基于tensorflow 概率计算。

There have been many questions and uncertainty around the future of PyMC3 since Theano stopped getting developed by the original authors, and we started experiments with a PyMC version based on tensorflow probability. 自从原作者停止开发 Theano 以来，PyMC3 的未来存在许多问题和不确定性，我们开始使用基于张量流概率的 PyMC 版本进行实验。

来源：pymc3·PyPI — pymc3 · PyPI

pymc与pymc3的安装与使用_pymc和pymc3-CSDN博客的帖子是2020.11的，帖子中选择安装pymc3可能是当时pymc还没迁移好，pymc3的Theano不能用可能也是历史原因。总之我认为现在应该安装pymc。

安装

官方推荐：安装 — PyMC 开发文档 — Installation — PyMC dev documentation

官方只给了用conda创建包含pymc的python环境的示例，我尝试在3.9python执行pip命令，发现也可以成功安装。

pip install pymc

PyMC中的变分推断

PyMC 支持各种变分推断技术。虽然这些方法速度更快，但它们通常也不太准确，并且可能导致推断有偏差。主要入口点是pymc.fit() 。

with pm.Model() as model:
    mu = pm.Normal("mu", mu=0, sigma=1)
    sd = pm.HalfNormal("sd", sigma=1)
    obs = pm.Normal("obs", mu=mu, sigma=sd, observed=rng.standard_normal(100))

    approx = pm.fit()

返回的Approximation对象具有各种功能，例如从近似后验中提取样本，我们可以像常规采样运行一样对其进行分析：

idata = approx.sample(1000)
az.summary(idata)

variational分子模块为 VI 的使用提供了很大的灵活性，并遵循面向对象的设计。例如，满秩 ADVI 估计完整协方差矩阵：

mu = pm.floatX([0.0, 0.0])
cov = pm.floatX([[1, 0.5], [0.5, 1.0]])
with pm.Model(coords={"idx": np.arange(2)}) as model:
    pm.MvNormal("x", mu=mu, cov=cov, dims="idx")
    approx = pm.fit(method="fullrank_advi")
'''上下等效'''
with pm.Model(coords={"idx": np.arange(2)}) as model:
    pm.MvNormal("x", mu=mu, cov=cov, dims="idx")
    approx = pm.FullRankADVI().fit()

Stein 变分梯度下降 (SVGD) 使用粒子来估计后验：

w = pm.floatX([0.2, 0.8])
mu = pm.floatX([-0.3, 0.5])
sd = pm.floatX([0.1, 0.1])
with pm.Model() as model:
    pm.NormalMixture("x", w=w, mu=mu, sigma=sd)
    approx = pm.fit(method=pm.SVGD(n_particles=200, jitter=1.0))
'''上下等效'''
with pm.Model() as model:
    pm.NormalMixture("x", w=[0.2, 0.8], mu=[-0.3, 0.5], sigma=[0.1, 0.1])

Pathfinder 变分推断

来源：Pathfinder 变分推理 — PyMC 示例库 — Pathfinder Variational Inference — PyMC example gallery

Pathfinder [张等人。 , 2021 ]是一种变分推理算法，可从贝叶斯模型的后验生成样本。它与广泛使用的 ADVI 算法相比毫不逊色。在大型问题上，它应该比大多数 MCMC 算法（包括动态 HMC（即 NUTS））具有更好的扩展性，但代价是对后验估计的偏差更大。有关算法的详细信息，请参阅arxiv 预印本。

该算法在BlackJAX （ JAX推理算法库）中实现。通过 PyMC 的 JAX 后端（通过pytensor ），我们可以使用一些简单的包装器代码在任何 PyMC 模型上运行 BlackJAX 的探路者。

此包装器代码在pymc-experimental中实现。本教程展示如何在 PyMC 模型上运行 Pathfinder。

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pymc as pm
import pymc_experimental as pmx

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

# Data of the Eight Schools Model
J = 8
y = np.array([28.0, 8.0, -3.0, 7.0, -1.0, 1.0, 18.0, 12.0])
sigma = np.array([15.0, 10.0, 16.0, 11.0, 9.0, 11.0, 10.0, 18.0])

with pm.Model() as model:
    mu = pm.Normal("mu", mu=0.0, sigma=10.0)
    tau = pm.HalfCauchy("tau", 5.0)

    z = pm.Normal("z", mu=0, sigma=1, shape=J)
    theta = mu + tau * z
    obs = pm.Normal("obs", mu=theta, sigma=sigma, shape=J, observed=y)

接下来，我们调用pmx.fit()并传入我们希望它使用的算法。

with model:
    idata = pmx.fit(method="pathfinder", num_samples=1000)

就像pymc.sample()一样，这会返回一个包含后验样本的 idata。请注意，由于这些样本并非来自 MCMC 链，因此无法以常规方式评估收敛性。

az.plot_trace(idata)
plt.tight_layout();

MLDA 采样器

来源：MLDA 采样器 — PyMC 示例库 — The MLDA sampler — PyMC example gallery

MLDA 采样器旨在处理计算密集型问题，我们不仅可以访问所需的（精细）后验分布，还可以访问一组降低精度和计算成本的近似（粗略）后验分布（我们至少需要以下之一）那些）。其主要思想是使用较粗链的样本作为较细链的建议。粗链运行固定次数的迭代，最后一个样本用作更细链的建议。这已被证明可以提高最细链的有效样本量，这使我们能够减少昂贵的细链似然评估的数量。

MLDA 的使用方式与 PyMC 中大多数步骤方法类似。它有一个特殊要求，即用户需要提供至少一个粗略模型才能使其工作。

'''定义精细模型fine_model'''
# Constructing the fine model
with pm.Model() as fine_model:
    # Define priors
    intercept = pm.Normal("intercept", 0, sigma=20)
    slope = pm.Normal("slope", 0, sigma=20)

    # Define likelihood
    likelihood = pm.Normal("y", mu=intercept + slope * x, sigma=sigma, observed=y)

    
'''定义粗略模型coarse_model'''
# Constructing the coarse model
with pm.Model() as coarse_model:
    # Define priors
    intercept = pm.Normal("intercept", 0, sigma=20)
    slope = pm.Normal("slope", 0, sigma=20)

    # Define likelihood
    likelihood = pm.Normal("y", mu=intercept + slope * x_coarse, sigma=sigma, observed=y_coarse)
  

'''使用 MLDA 从后验中绘制 MCMC 样本'''
with fine_model:
    # Initialise step methods
    step = pm.MLDA(coarse_models=[coarse_model], subsampling_rates=[10])
    step_2 = pm.Metropolis()
    step_3 = pm.DEMetropolisZ()

    # Sample using MLDA
    t_start = time.time()
    trace = pm.sample(draws=6000, chains=4, tune=2000, step=step, random_seed=RANDOM_SEED)
    runtime = time.time() - t_start

    # Sample using Metropolis
    t_start = time.time()
    trace_2 = pm.sample(draws=6000, chains=4, tune=2000, step=step_2, random_seed=RANDOM_SEED)
    runtime_2 = time.time() - t_start

    # Sample using DEMetropolisZ
    t_start = time.time()
    trace_3 = pm.sample(draws=6000, chains=4, tune=2000, step=step_3, random_seed=RANDOM_SEED)
    runtime_3 = time.time() - t_start

我们将coarse_model提供给 MLDA 实例，并将subsampling_rate设置为 10。子采样率是在粗链中抽取的样本数量，用于构建细链的提案。在这种情况下，MLDA 在粗链中抽取 10 个样本，并使用最后一个作为细链的提案。这被细链接受或拒绝，然后控制返回到粗链，粗链生成另外 10 个样本等。请注意， pm.MLDA还有许多其他调整参数，可以在文档中找到。

接下来，我们使用通用pm.sample方法，将 MLDA 实例传递给它。这会运行 MLDA 并返回trace ，其中包含所有 MCMC 样本和各种副产品。在这里，我们还运行标准 Metropolis 和 DEMetropolisZ 采样器进行比较，它们返回单独的跟踪。我们对运行进行计时以便稍后进行比较。

最后，PyMC 提供了各种函数来可视化跟踪和打印摘要统计信息（其中两个如下所示）。

# Trace plots
az.plot_trace(trace)
az.plot_trace(trace_2)
az.plot_trace(trace_3)

经验近似概述

对于大多数模型，我们使用采样 MCMC 算法，例如 Metropolis 或 NUTS。在 PyMC 中，我们习惯于存储 MCMC 样本的痕迹，然后使用它们进行分析。 PyMC 中的变分推理子模块有一个类似的概念： Empirical 。这种类型的近似存储 SVGD 采样器的粒子。独立的SVGD粒子和MCMC样本之间没有区别。 Empirical充当 MCMC 采样输出和成熟的 VI 实用程序（如apply_replacements或sample_node之间的桥梁。接口说明请参见variational_api_quickstart 。在这里，我们将只关注Emprical并概述经验近似的具体内容。

import arviz as az
import matplotlib.pyplot as plt
import numpy as np
import pymc as pm
import pytensor
import seaborn as sns

from pandas import DataFrame

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

对于variational_api_quickstart中的问题，首先使用默认的NUTS求解：

w = pm.floatX([0.2, 0.8])
mu = pm.floatX([-0.3, 0.5])
sd = pm.floatX([0.1, 0.1])

with pm.Model() as model:
    x = pm.NormalMixture("x", w=w, mu=mu, sigma=sd)
    trace = pm.sample(50_000, return_inferencedata=False)
    
with model:
    idata = pm.to_inference_data(trace)
az.plot_trace(idata);

如果要创建Empirical近似值：

with model:
    approx = pm.Empirical(trace)

这种类型的近似有自己的样本底层存储，即pytensor.shared本身。它的样本数量与您之前跟踪的样本数量完全相同。在我们的当前情况下，它是 50k。另一件需要注意的事情是，如果您的多踪迹具有多个链，您将同时存储更多的样本。我们将所有痕迹展平以创建Empirical 。

后验采样是通过替换统一完成的。调用approx.sample(1000) ，您将再次获得跟踪，但顺序尚未确定。但是之后无法使用approx.sample再次重建底层跟踪。

编译采样函数后，采样速度变得非常快。您会看到不再有顺序，但重建的密度是相同的。

new_trace = approx.sample(50000)
az.plot_trace(new_trace);

结语

下期预告：变分推断：贝叶斯神经网络

概率编程、深度学习和“大数据”是机器学习中最热门的话题。在概率编程中，很多创新都集中在使用变分推断来扩展事物。在此示例中，我将展示如何在 PyMC 中使用变分推断来拟合简单的贝叶斯神经网络。我还将讨论如何将概率编程和深度学习结合起来，为未来的研究开辟非常有趣的探索途径。

[[外链图片转存中...(img-muN8POk9-1725243832234)]](https://postimg.cc/SjRFtNDw)

## 结语

下期预告：变分推断：贝叶斯神经网络

**概率编程**、**深度学习**和“**大数据**”是机器学习中最热门的话题。在概率编程中，很多创新都集中在使用**变分**推断来扩展事物。在此示例中，我将展示如何在 PyMC 中使用**变分推断**来拟合简单的贝叶斯神经网络。我还将讨论如何将概率编程和深度学习结合起来，为未来的研究开辟非常有趣的探索途径。

欢迎在评论区交流讨论，某乎、某书、某号同名，期待您的关注！