​Pyro简介 贝叶斯神经网络bnn , 隐马尔可夫模型 人工智能python python 概率分布程序包的使用教程

Pyro简介

镜像GitCode - 全球开发者的开源社区,开源代码托管平台

gitcode.com/gh_mirrors/py/pyro/blob/dev/tutorial/source/bayesian_regression_ii.ipynb

Introduction to Pyro — Pyro Tutorials 1.9.1 documentation

pyro.ai/examples/intro_long.html#Example-model:-Maximum-likelihood-linear-regression

概率是在不确定性下推理的数学，就像微积分是推理变化率的数学一样。它为理解许多现代机器学习和人工智能提供了一个统一的理论框架:用概率语言建立的模型可以捕捉复杂的推理，知道他们不知道的东西，并在没有监督的情况下揭示数据的结构。

直接指定概率模型可能很麻烦，并且实现它们可能很容易出错。概率编程语言(ppl)通过将概率与编程语言的表达能力相结合来解决这些问题。概率程序是普通确定性计算和随机采样值的混合，表示生成过程为了数据。

通过观察概率程序的结果，我们可以描述一个推理问题，大致翻译为:“如果这个随机选择有一定的观察值，那么什么一定是真的？”ppl明确地强制实施了一种关注点的分离，这种分离已经隐含在模型规范、要回答的查询和计算答案的算法之间的概率数学中。

Pyro是一种基于Python和PyTorch的概率编程语言。Pyro程序只是Python程序，而它的主要推理技术是随机变分推理，它将抽象的概率计算转化为用PyTorch中的随机梯度下降解决的具体优化问题，使概率方法适用于以前难以处理的模型和数据集大小。

在本教程中，我们对概率机器学习和用Pyro进行概率编程的基本概念进行了一个简短的、自以为是的浏览。我们通过一个涉及线性回归的示例数据分析问题来做到这一点，线性回归是机器学习中最常见和最基本的任务之一。我们将看到如何使用Pyro的建模语言和推理算法将不确定性合并到回归系数的估计中。

概述¶

• 介绍

• 概述

• 设置

• 背景:概率机器学习

  • 背景:概率模型

  • 背景:推理、学习和评估

• 例如:地理和国民收入

• Pyro中的模型

  • 示例模型:最大似然线性回归

  • 背景:pyro.sample原语

  • 背景:pyro.param原语

  • 背景:pyro.plate原语

  • 示例:从最大似然回归到贝叶斯回归

• 烟火中的推理

  • 背景:变分推理

  • 背景:“指南”程序作为灵活的近似后验概率

  • 示例:Pyro中贝叶斯线性回归的平均场变分近似

  • 背景:估计和优化证据下限(ELBO)

  • 示例:通过随机变分推理的贝叶斯回归(SVI)

• Pyro中的模型评估

  • 背景:贝叶斯模型评估与后验预测检查

  • 示例:Pyro中的后验预测不确定性

  • 示例:用全秩指南重新审视贝叶斯回归

• 后续步骤

设置¶

让我们从导入我们需要的模块开始。

[41]:

%reset -s -f

[42]:

import logging
import os

import torch
import numpy as np
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt

import pyro

[43]:

smoke_test = ('CI' in os.environ)
assert pyro.__version__.startswith('1.9.1')

pyro.enable_validation(True)
pyro.set_rng_seed(1)
logging.basicConfig(format='%(message)s', level=logging.INFO)

# Set matplotlib settings
%matplotlib inline
plt.style.use('default')

背景:概率机器学习¶

大多数数据分析问题可以理解为对三个基本高层次问题的阐述:

1. 在观察任何数据之前，我们对问题有什么了解？

2. 根据我们的先验知识，我们可以从数据中得出什么结论？

3. 这些结论有意义吗？

在数据科学和机器学习的概率或贝叶斯方法中，我们根据概率分布的数学运算来形式化这些。

背景:概率模型¶

首先，我们把我们所知道的关于问题中变量的一切以及它们之间的关系用概率模型或者随机变量集合上的联合概率分布。一个模型有观察x和潜在的随机变量z以及参数θ。它通常具有以下形式的联合密度函数

pθ(x,z)=pθ(x|z)pθ(z)

潜在变量的分布pθ(z)在这个公式中称为在先的；在前的，以及给定潜在变量的观察变量的分布pθ(x|z)被称为可能性.

我们通常要求各种条件概率分布pi组成了一个模型pθ(x,z)具有以下属性(通常满足中提供的分布)放火狂者和PyTorch分布):

• 我们可以有效地从每个样本中pi

• 我们可以有效地计算逐点概率密度pi

• pi关于参数是可微的θ

概率模型通常用标准图形符号用于可视化和交流，总结如下，尽管在Pyro中可以表示没有固定图形结构的模型。在有很多重复的模型中，使用起来很方便盘子符号，之所以这样叫是因为它以图形方式显示为一个围绕变量的矩形“盘子”,以表示里面随机变量的多个独立副本。

背景:推理、学习和评估¶

一旦我们指定了一个模型，贝叶斯法则告诉我们如何使用它来执行推理，或者从数据中得出关于潜在变量的结论后验分布超过z:

pθ(z|x)=pθ(x,z)∫dzpθ(x,z)

为了检查建模和推理的结果，我们想知道模型与观察到的数据有多吻合x，我们可以用证据或者边际可能性

pθ(x)=∫dzpθ(x,z)

还可以预测新的数据，我们可以用后验预测分布

pθ(x′|x)=∫dzpθ(x′|z)pθ(z|x)

最后，通常希望学习参数θ我们的模型来自观测数据x，我们可以通过最大化边际可能性来实现:

θmax=argmaxθpθ(x)=argmaxθ∫dzpθ(x,z)

例如:地理和国民收入¶

下面的例子改编自这本优秀的书的第七章统计学反思作者Richard McElreath，鼓励读者阅读这篇文章，以获得对贝叶斯数据分析的更广泛实践的简单介绍(所有章节的火焰代码可用)。

我们想探索一个国家的地形异质性之间的关系，由地形粗糙度指数(变量崎岖的数据集中)及其人均GDP。特别是，原始论文的作者指出(“崎岖:非洲恶劣地理环境的祝福”地形崎岖或恶劣的地理环境与非洲以外的经济表现较差有关，但崎岖的地形对非洲国家的收入产生了相反的影响。

让我们看看数据，研究一下这种关系。我们将重点关注数据集中的三个要素:

• rugged:量化地形崎岖指数；

• cont_africa给定民族是否在非洲；

• rgdppc_2000:2000年实际人均国内生产总值；

[44]:

DATA_URL = "https://d2hg8soec8ck9v.cloudfront.net/datasets/rugged_data.csv"
data = pd.read_csv(DATA_URL, encoding="ISO-8859-1")
df = data[["cont_africa", "rugged", "rgdppc_2000"]]

响应变量GDP是高度倾斜的，因此我们将在继续之前对其进行对数转换。

[45]:

df = df[np.isfinite(df.rgdppc_2000)]
df["rgdppc_2000"] = np.log(df["rgdppc_2000"])

然后，我们将这个数据帧后面的Numpy数组转换为torch.Tensor用于PyTorch和Pyro的分析。

[46]:

train = torch.tensor(df.values, dtype=torch.float)
is_cont_africa, ruggedness, log_gdp = train[:, 0], train[:, 1], train[:, 2]

可视化数据表明，坚固性和GDP之间确实存在可能的关系，但需要进一步的分析来证实这一点。我们将看到如何通过贝叶斯线性回归在Pyro中做到这一点。

[47]:

fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(12, 6), sharey=True)
african_nations = df[df["cont_africa"] == 1]
non_african_nations = df[df["cont_africa"] == 0]
sns.scatterplot(x=non_african_nations["rugged"],
                y=non_african_nations["rgdppc_2000"],
                ax=ax[0])
ax[0].set(xlabel="Terrain Ruggedness Index",
          ylabel="log GDP (2000)",
          title="Non African Nations")
sns.scatterplot(x=african_nations["rugged"],
                y=african_nations["rgdppc_2000"],
                ax=ax[1])
ax[1].set(xlabel="Terrain Ruggedness Index",
          ylabel="log GDP (2000)",
          title="African Nations");

Pyro中的模型¶

Pyro中的概率模型被指定为Python函数model(*args, **kwargs)从潜在变量中产生观察数据特殊原函数它的行为可以由Pyro的内部根据正在执行的高级计算来改变。

具体来说，不同的数学部分model()通过映射进行编码:

1. 潜在随机变量z ⟺ 高温样品

2. 观察随机变量x ⟺ 高温样品与obs关键字参数

3. 可学习的参数θ ⟺ pyro.param

4. 盘子⟺ 焦碳板上下文管理器

我们在下面的线性回归的第一个Pyro模型的上下文中详细检查每个组件。

示例模型:最大似然线性回归¶

如果我们写出线性回归预测值的公式βX+α作为Python表达式，我们得到以下内容:

mean = a + b_a * is_cont_africa + b_r * ruggedness + b_ar * is_cont_africa * ruggedness

为了将它构建成数据集的完整概率模型，我们需要确定回归系数α和β可学习的参数并在预测平均值周围添加观察噪声。我们可以使用上面介绍的Pyro原语来表达这一点，并使用pyro.render_model():

[48]:

import pyro.distributions as dist
import pyro.distributions.constraints as constraints

def simple_model(is_cont_africa, ruggedness, log_gdp=None):
    a = pyro.param("a", lambda: torch.randn(()))
    b_a = pyro.param("bA", lambda: torch.randn(()))
    b_r = pyro.param("bR", lambda: torch.randn(()))
    b_ar = pyro.param("bAR", lambda: torch.randn(()))
    sigma = pyro.param("sigma", lambda: torch.ones(()), constraint=constraints.positive)

    mean = a + b_a * is_cont_africa + b_r * ruggedness + b_ar * is_cont_africa * ruggedness

    with pyro.plate("data", len(ruggedness)):
        return pyro.sample("obs", dist.Normal(mean, sigma), obs=log_gdp)

pyro.render_model(simple_model, model_args=(is_cont_africa, ruggedness, log_gdp), render_distributions=True)

[48]:

上面的图没有显示模型参数a, b_a, b_r, b_ar，以及sigma。我们可以设定render_params=True渲染模型参数。

[49]:

pyro.render_model(simple_model, model_args=(is_cont_africa, ruggedness, log_gdp), render_distributions=True, render_params=True)

[49]:

论点render_distributions = True将显示对参数的约束。举个例子，sigma是非负的标准差。因此，其约束显示为“sigma∈GreaterThan(lower_bound=0.0)"。

学习的参数simple_model会构成最大似然估计并产生回归系数的点估计。然而，在这个例子中，我们的数据可视化表明，我们不应该对回归系数的任何单一值过于自信。相比之下，完全贝叶斯方法将对不同的可能参数值以及模型预测产生不确定性估计。

在制作线性模型的贝叶斯版本之前，让我们暂停一下，仔细看看第一段Pyro代码。

背景pyro.sample原始的¶

Pyro中的概率程序是围绕原始概率分布的样本构建的，标记为pyro.sample:

def sample(
    name: str,
    fn: pyro.distributions.Distribution,
    *,
    obs: typing.Optional[torch.Tensor] = None,
    infer: typing.Optional[dict] = None
) -> torch.Tensor:
    ...

在我们的模型中simple_model在这条线上

return pyro.sample("obs", dist.Normal(mean, sigma), obs=log_gdp)

可以代表潜在变量或观察变量，取决于是否simple_model被赋予一个log_gdp价值。当...的时候log_gdp未提供，并且"obs"是潜在的，它相当于(忽略了pyro.plate现在通过假设len(ruggedness) == 1)来调用发行版的底层.sample方法:

return dist.Normal(mean, sigma).sample()

这种解释是偶尔提到烟火计划的原因随机函数在Pyro的一些旧文档中使用了一个相当模糊的术语。

当...的时候simple_model被赋予一个log_gdp参数和"obs"是观察到的pyro.sample声明

return pyro.sample("obs", dist.Normal(mean, sigma), obs=log_gdp)

总会回来的log_gdp:

return log_gdp

但是，请注意当观察到任何示例语句时，模型中每个其他示例语句的累积效果都会发生变化遵循贝叶斯法则；这是派若的工作推理算法“逆向运行程序”并为所有程序分配数学上一致的值pyro.sample模型中的语句。

此时，一个合理的问题是为什么pyro.sample而其他原语必须有名字。用户和Pyro的内部利用名称来分离模型、观察和推理算法的规范，这是概率编程语言的一个关键卖点。为了看到这样的例子，我们可以看看高阶原语高温条件这解决了针对单个Pyro模型编写许多查询的问题:

def condition(
    model: Callable[..., T],
    data: Dict[str, torch.Tensor]
) -> Callable[..., T]:
    ...

pyro.condition获取一个模型和一个(可能为空的)字典，将名称映射到观察值，并将每个观察值传递给pyro.sample由其名称表示的语句。在我们的例子的上下文中simple_model我们可以移除log_gdp作为一个参数，并用一个更简单的接口代替它

def simpler_model(is_cont_africa, ruggedness): ...

conditioned_model = pyro.condition(simpler_model, data={"obs": log_gdp})

在哪里conditioned_model相当于

conditioned_model = functools.partial(simple_model, log_gdp=log_gdp)

背景pyro.param原始的¶

我们模型中使用的下一个原语，pyro.param是读取和写入Pyro的前端键值参数存储:

def param(
    name: str,
    init: Optional[Union[torch.Tensor, Callable[..., torch.Tensor]]] = None,
    *,
    constraint: torch.distributions.constraints.Constraint = constraints.real
) -> torch.Tensor:
    ...

喜欢pyro.sample, pyro.param总是以名称作为第一个参数来调用。第一次pyro.param用特定的名称调用时，它存储由第二个参数指定的初始值init然后返回该值。此后，当使用该名称调用它时，它将从参数存储中返回该值，而不考虑任何其他参数。参数初始化后，不再需要指定init为了检索其值(例如pyro.param("a")).

第二个论点，init，可以是torch.Tensor或者一个不带参数但返回张量的函数。第二种形式很有用，因为它避免了重复构造仅在模型第一次运行时使用的初始值。

不像PyTorch的torch.nn.Parameter的情况下，Pyro中的参数可以显式地约束到Rn这是一个重要的特征，因为许多基本概率分布的参数具有受限的域。例如，在scale的参数Normal分配必须是正的。的可选第三个参数pyro.param, constraint，是一个火炬.分布.约束.约束初始化参数时存储的对象；每次更新后都会重新应用约束。Pyro附带了大量预定义的约束。

pyro.param除非参数存储区由优化算法更新或通过pyro.clear_param_store()。不像pyro.sample, pyro.param在一个模型中可以用相同的名称调用多次；每个同名的调用都将返回相同的值。全局参数存储本身可通过调用pyro.get_param_store().

在我们的模型中simple_model以上，声明

a = pyro.param("a", lambda: torch.randn(()))

在概念上类似于下面的代码，但是增加了一些跟踪、序列化和约束管理功能。

simple_param_store = {}
...
def simple_model():
    a_init = lambda: torch.randn(())
    a = simple_param_store["a"] if "a" in simple_param_store else a_init()
    ...

虽然本介绍性教程使用了pyro.param对于参数管理，Pyro也兼容PyTorch的familiartorch.nn.ModuleAPI via高温模块.

背景pyro.plate原始的¶

焦碳板是Pyro的正式编码平板符号，广泛用于概率机器学习中，以简化模型的可视化和分析独立同分布随机变量。

def plate(
    name: str,
    size: int,
    *,
    dim: Optional[int] = None,
    **other_kwargs
) -> contextlib.AbstractContextManager:
    ...

概念上，pyro.plate语句相当于for-循环。在…里simple_model我们可以替换这些线条

with pyro.plate("data", len(ruggedness)):
    return pyro.sample("obs", dist.Normal(mean, sigma), obs=log_gdp)

用一条蟒蛇for-循环

result = torch.empty_like(ruggedness)
for i in range(len(ruggedness)):
    result[i] = pyro.sample(f"obs_{i}", dist.Normal(mean, sigma), obs=log_gdp[i] if log_gdp is not None else None)
return result

当重复变量的数量(len(ruggedness)在这种情况下)很大，使用Python循环会非常慢。因为循环中的每次迭代都是相互独立的，pyro.plate使用PyTorch的阵列广播要在单个矢量化操作中并行执行迭代，如下面的等价矢量化代码所示:

mean = mean.unsqueeze(-1).expand((len(ruggedness),))
sigma = sigma.unsqueeze(-1).expand((len(ruggedness),))
return pyro.sample("obs", dist.Normal(mean, sigma), obs=log_gdp)

实际上，在编写Pyro程序时pyro.plate是最有用的矢量化工具。但是，如中所述SVI教程的第2部分，它也有助于管理数据子采样并且作为推断算法的信号，表明某些变量是独立的。

示例:从最大似然回归到贝叶斯回归¶

为了使我们的线性回归模型贝叶斯化，我们需要指定先验分布关于参数α∈R和β∈R3(此处扩展为标量b_a, b_r，以及b_ar).这些概率分布代表了我们在观察任何关于合理值的数据之前的信念α和β。我们还将添加一个随机比例参数σ控制观察噪音。

在Pyro中表达线性回归的贝叶斯模型非常直观:我们只需替换每个pyro.param语句与pyro.sample报表配有烟火分配物体描述每个参数的先验信念。

对于常数项α，我们使用一个标准偏差较大的正态先验来表示我们相对缺乏关于基线GDP的先验知识。对于其他回归系数，我们使用标准正态先验(以0为中心)来表示我们缺乏关于协变量和GDP之间的关系是正还是负的先验知识。对于观察噪声σ我们使用一个低于0的平坦先验，因为这个值必须是正的才是有效的标准偏差。

[50]:

def model(is_cont_africa, ruggedness, log_gdp=None):
    a = pyro.sample("a", dist.Normal(0., 10.))
    b_a = pyro.sample("bA", dist.Normal(0., 1.))
    b_r = pyro.sample("bR", dist.Normal(0., 1.))
    b_ar = pyro.sample("bAR", dist.Normal(0., 1.))
    sigma = pyro.sample("sigma", dist.Uniform(0., 10.))

    mean = a + b_a * is_cont_africa + b_r * ruggedness + b_ar * is_cont_africa * ruggedness

    with pyro.plate("data", len(ruggedness)):
        return pyro.sample("obs", dist.Normal(mean, sigma), obs=log_gdp)

pyro.render_model(model, model_args=(is_cont_africa, ruggedness, log_gdp), render_distributions=True)

[50]:

烟火中的推理¶

背景:变分推理¶

简介中的每个计算(后验分布、边际可能性和后验预测分布)都需要进行积分，而这些积分通常是不可能的或计算上难以处理的。

虽然Pyro支持许多不同的精确和近似推理算法，但最受支持的是变分推理，它为查找提供了一个统一的方案θmax以及计算易处理的近似值qϕ(z)到真实的，未知的后路pθmax(z|x)通过将难以处理的积分转化为p和q。下图从概念上描述了这一过程，而更全面的数学介绍可在SVI教程.

大多数概率分布(下图中的浅椭圆)，尤其是对应于贝叶斯后验分布的概率分布，太复杂而无法直接表示，因此我们必须定义一个更小的子空间，由实值参数索引ϕ分配的qϕ(z)从结构上保证易于取样(下图中的黑圈)，但可能不包括真实的后验分布pθ(z|x)(下图中的红星)。

变分推理通过搜索变分分布的空间来近似真实的后验概率，以根据一些距离或散度的度量(下图中的黑色箭头)找到与真实的后验概率(下图中的黄色星)最相似的一个。

然而，有许多不同的方法来测量概率分布之间的距离或散度。我们应该选择哪一个？如图所示，理论上有吸引力的选择是库尔贝克-莱布勒散度KL(qϕ(z)||pθ(z|x))，但是直接计算这需要提前知道真实的后验概率，这将违背目的。

更重要的是，我们对优化这种分歧感兴趣，这听起来可能更难，但事实上，使用贝叶斯定理来重写的定义是可能的KL(qϕ(z)||pθ(z|x))一个棘手的常数和一个不依赖于qϕ一个容易理解的术语叫做证据下限(ELBO)，定义如下。因此，最大化这个易处理项将产生与最小化原始KL散度相同的解。

背景:“指南”程序作为灵活的近似后验概率¶

在变分推断中，我们引入了一个参数化分布qϕ(z)接近真实的后部，在哪里ϕ被称为变分参数。这个分布在很多文献中被称为变分分布，在Pyro的上下文中被称为向导(一个音节而不是九个！).

就像模型一样，指南被编码为Python程序guide()其中包含pyro.sample和pyro.param声明。确实如此不包含观察到的数据，因为指南需要是一个适当的标准化分布，以便于取样。请注意，Pyro强制执行这一点model()和guide()应该采取同样的观点。允许向导是任意的Pyro程序打开了编写向导族的可能性，这些向导族捕获更多的真实后验概率的特定问题结构，只在有帮助的方向上扩展搜索空间，如下图所示。

变分推理的数学对指南施加了什么限制？因为导向器是后部的近似值pθmax(z|x)，指南需要提供模型中所有潜在随机变量的有效联合概率密度。回想一下，当在Pyro中用原语语句指定随机变量时pyro.sample()第一个参数表示随机变量的名称。这些名称将用于对齐模型和指南中的随机变量。非常明确地说，如果模型包含一个随机变量z_1

def model():
    pyro.sample("z_1", ...)

那么向导需要有一个匹配sample声明

def guide():
    pyro.sample("z_1", ...)

这两种情况下使用的分布可以不同，但是名称必须一一对应。

尽管它提供了很大的灵活性，但是手工编写指南可能会很困难和乏味，尤其是对于新用户来说。只要有可能，我们建议使用自动制导，或用于从附带Pyro的模型中自动生成公共导轨族的配方自动制导。下一节将演示这两种方法。

示例:Pyro中贝叶斯线性回归的平均场变分近似¶

对于贝叶斯线性回归的运行示例，我们将使用一个指南，将模型中未观察到的参数的分布建模为具有对角协方差的高斯分布，即假设潜在变量之间没有相关性(我们将看到这是一个非常强的假设)。这被称为平均场近似这是一个从物理学借用的术语，这种近似法最初就是在物理学中发明的。

为了完整起见，我们先用手写出这种形式的引导程序。

[51]:

def custom_guide(is_cont_africa, ruggedness, log_gdp=None):
    a_loc = pyro.param('a_loc', lambda: torch.tensor(0.))
    a_scale = pyro.param('a_scale', lambda: torch.tensor(1.),
                         constraint=constraints.positive)
    sigma_loc = pyro.param('sigma_loc', lambda: torch.tensor(0.))
    weights_loc = pyro.param('weights_loc', lambda: torch.randn(3))
    weights_scale = pyro.param('weights_scale', lambda: torch.ones(3),
                               constraint=constraints.positive)
    a = pyro.sample("a", dist.Normal(a_loc, a_scale))
    b_a = pyro.sample("bA", dist.Normal(weights_loc[0], weights_scale[0]))
    b_r = pyro.sample("bR", dist.Normal(weights_loc[1], weights_scale[1]))
    b_ar = pyro.sample("bAR", dist.Normal(weights_loc[2], weights_scale[2]))
    sigma = pyro.sample("sigma", dist.LogNormal(sigma_loc, torch.tensor(0.05)))  # fixed scale for simplicity
    return {"a": a, "b_a": b_a, "b_r": b_r, "b_ar": b_ar, "sigma": sigma}

我们可以使用pyro.render_model想象custom_guide，证明了随机变量确实是相互独立的，因为它们之间没有边。

[52]:

pyro.render_model(custom_guide, model_args=(is_cont_africa, ruggedness, log_gdp), render_params=True)

[52]:

Pyro还包含一个广泛的“自动向导”集合，它可以根据给定的模型自动生成向导程序。就像我们手写的指南一样pyro.autoguide.AutoGuide实例(本身就是接受与模型相同的参数的函数)返回每个值的字典pyro.sample它们包含的站点。

最简单的自动引导类是AutoNormal，它在一行代码中自动生成一个指南，相当于我们上面手写的指南:

[53]:

auto_guide = pyro.infer.autoguide.AutoNormal(model)

然而，指南本身并没有完全指定推理算法:它仅仅描述了由参数(上图中的黑圈)索引的可能近似后验分布上的搜索空间，以及由初始参数值确定的该空间中的初始点。然后，我们必须通过解决参数的优化问题(上图中的黄色星号)，将初始分布移向真实的后验分布(上图中的红色星号)。制定和解决这个优化问题是下两节的主题。

背景:估计和优化证据下限(ELBO)¶

模型的功能pθ(x,z)和引导qϕ(z)我们将优化的是ELBO，它被定义为指南中样本的预期w.r.t:。

ELBO≡Eqϕ(z)[原木⁡pθ(x,z)−原木⁡qϕ(z)]

通过假设，我们可以计算期望中的所有概率，因为指南q假设是我们可以从中采样的参数分布，我们可以计算这个量以及关于模型和引导参数的梯度的蒙特卡罗估计，∇θ,ϕELBO.

优化ELBO over模型和指南参数θ,ϕ通过随机梯度下降使用这些梯度估计有时被称为随机变分推理(SVI)；有关SVI的详细介绍，请参阅SVI第一部分.

示例:通过随机变分推理的贝叶斯回归(SVI)¶

Pyro含有许多不同的实现方式ELBO的估计器(在前面的章节中进行了数学定义),每个估计器计算损耗和梯度的方式稍有不同。在本教程中，我们将只使用pyro.infer.Trace_ELBO，这样做总是正确和安全的；其他ELBO估计量可能为某些模型和指南提供计算或统计优势。

我们将在示例模型中使用SVI进行推理，演示Pyro如何使用PyTorch的随机梯度下降实现来优化pyro.infer.Trace_ELBO我们传递给的对象pyro.infer.SVI，这是一个帮助器类，其step()方法负责计算损失和参数梯度，并对参数应用更新和约束。

[54]:

adam = pyro.optim.Adam({"lr": 0.02})
elbo = pyro.infer.Trace_ELBO()
svi = pyro.infer.SVI(model, auto_guide, adam, elbo)

这里pyro.optim.Adam是PyTorch优化器的一层薄薄的包装火炬. optim .亚当(参见这里供讨论)。优化器在pyro.optim用于优化和更新Pyro参数存储中的参数值。特别是，您会注意到我们不需要向优化器传递可学习的参数，因为这是由指导代码决定的，并且在SVI自动分类。采取埃尔伯梯度步骤，我们简单地称之为SVI的步骤方法。我们传递给的数据参数SVI.step将传递给双方model()和guide()。完整的训练循环如下:

[55]:

%%time
pyro.clear_param_store()

# These should be reset each training loop.
auto_guide = pyro.infer.autoguide.AutoNormal(model)
adam = pyro.optim.Adam({"lr": 0.02})  # Consider decreasing learning rate.
elbo = pyro.infer.Trace_ELBO()
svi = pyro.infer.SVI(model, auto_guide, adam, elbo)

losses = []
for step in range(1000 if not smoke_test else 2):  # Consider running for more steps.
    loss = svi.step(is_cont_africa, ruggedness, log_gdp)
    losses.append(loss)
    if step % 100 == 0:
        logging.info("Elbo loss: {}".format(loss))

plt.figure(figsize=(5, 2))
plt.plot(losses)
plt.xlabel("SVI step")
plt.ylabel("ELBO loss");

Elbo loss: 694.9404826164246
Elbo loss: 524.3822101354599
Elbo loss: 475.66820669174194
Elbo loss: 399.99088364839554
Elbo loss: 315.23274326324463
Elbo loss: 254.76771265268326
Elbo loss: 248.237040579319
Elbo loss: 248.42670530080795
Elbo loss: 248.46450632810593
Elbo loss: 257.41463351249695

CPU times: user 6.47 s, sys: 241 µs, total: 6.47 s
Wall time: 6.28 s

[55]:

Text(0, 0.5, 'ELBO loss')

请注意，由于我们使用了较高的学习率，因此本次培训速度很快。有时模型和向导对学习速度很敏感，首先要尝试的是降低学习速度和增加步数。这在具有深度神经网络的模型和向导中尤其重要。我们建议从较低的学习速度开始，然后逐渐增加，避免学习速度过快，否则推理会出现偏差或导致错误。

训练好向导后，我们可以通过从Pyro的param存储中获取来检查优化的向导参数值。每个(loc, scale)下面打印的对参数化单个pyro.distributions.Normal指南中的分布对应于不同的未观察到的pyro.sample模型中的语句，类似于我们的手写custom_guide以前的。

[56]:

for name, value in pyro.get_param_store().items():
    print(name, pyro.param(name).data.cpu().numpy())

AutoNormal.locs.a 9.173145
AutoNormal.scales.a 0.0703669
AutoNormal.locs.bA -1.8474661
AutoNormal.scales.bA 0.1407009
AutoNormal.locs.bR -0.19032118
AutoNormal.scales.bR 0.044044234
AutoNormal.locs.bAR 0.35599768
AutoNormal.scales.bAR 0.079374395
AutoNormal.locs.sigma -2.205863
AutoNormal.scales.sigma 0.060526706

最后，让我们重温一下我们之前的问题，即地形崎岖度和GDP之间的关系对于模型参数估计中的任何不确定性有多强。为此，我们绘制了给定地形崎岖度的非洲内外国家的对数GDP的斜率分布。

我们用从我们训练有素的向导中抽取的样本来表示这两种分布。要并行绘制多个样本，我们可以在pyro.plate语句，该语句重复并向量化每个pyro.sample语句，如介绍pyro.plate原始。

[57]:

with pyro.plate("samples", 800, dim=-1):
    samples = auto_guide(is_cont_africa, ruggedness)

gamma_within_africa = samples["bR"] + samples["bAR"]
gamma_outside_africa = samples["bR"]

从下面可以看出，非洲国家的概率主要集中在正区域，其他国家则相反，这进一步证实了最初的假设。然而，非非洲国家的后验不确定性(橙色直方图的宽度)似乎大大低于非洲国家(蓝色直方图的宽度)，这是令人惊讶的，因为原始数据中似乎有类似的分布。我们将在下一节进一步研究这种差异。

[58]:

fig = plt.figure(figsize=(10, 6))
sns.histplot(gamma_within_africa.detach().cpu().numpy(), kde=True, stat="density", label="African nations")
sns.histplot(gamma_outside_africa.detach().cpu().numpy(), kde=True, stat="density", label="Non-African nations", color="orange")
fig.suptitle("Density of Slope : log(GDP) vs. Terrain Ruggedness");
plt.xlabel("Slope of regression line")
plt.legend()
plt.show()

Pyro中的模型评估¶

背景:贝叶斯模型评估与后验预测检查¶

为了评估我们是否可以信任我们的推断结果，我们将比较由我们的模型诱导的可能新数据的后验预测分布与现有的观察数据。计算这种分布通常是困难的，因为它依赖于知道真实的后验概率，但是我们可以很容易地使用从变分推断中获得的近似后验概率来近似它:

pθ(x′|x)=∫dzpθ(x′|z)pθ(z|x)≈∫dzpθ(x′|z)qϕ(z|x)

具体地说，为了从后验预测中抽取一个近似样本，我们简单地抽取一个样本z^∼qϕ(z)从近似后验概率，然后从我们模型中给定样本的观察变量的分布中抽取样本x′∼pθ(x|z^)，好像我们用(近似的)后验代替了先验。

示例:Pyro中的后验预测不确定性¶

为了评估我们的示例线性回归模型，我们将使用预言性的实用程序类，它实现了上面的方法，大约从pθ(x′|x).

我们从训练好的模型中生成800个样本。在内部，这是通过首先从guide，然后向前运行模型，同时更改unobserved返回的值pyro.sample语句中采样的相应值guide.

[59]:

predictive = pyro.infer.Predictive(model, guide=auto_guide, num_samples=800)
svi_samples = predictive(is_cont_africa, ruggedness, log_gdp=None)
svi_gdp = svi_samples["obs"]

下面的代码是特定于这个例子的，只是用来绘制每个国家后验预测分布的90%可信区间(包含90%概率质量的区间)。

[60]:

predictions = pd.DataFrame({
    "cont_africa": is_cont_africa,
    "rugged": ruggedness,
    "y_mean": svi_gdp.mean(0).detach().cpu().numpy(),
    "y_perc_5": svi_gdp.kthvalue(int(len(svi_gdp) * 0.05), dim=0)[0].detach().cpu().numpy(),
    "y_perc_95": svi_gdp.kthvalue(int(len(svi_gdp) * 0.95), dim=0)[0].detach().cpu().numpy(),
    "true_gdp": log_gdp,
})
african_nations = predictions[predictions["cont_africa"] == 1].sort_values(by=["rugged"])
non_african_nations = predictions[predictions["cont_africa"] == 0].sort_values(by=["rugged"])

fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(12, 6), sharey=True)
fig.suptitle("Posterior predictive distribution with 90% CI", fontsize=16)

ax[0].plot(non_african_nations["rugged"], non_african_nations["y_mean"])
ax[0].fill_between(non_african_nations["rugged"], non_african_nations["y_perc_5"], non_african_nations["y_perc_95"], alpha=0.5)
ax[0].plot(non_african_nations["rugged"], non_african_nations["true_gdp"], "o")
ax[0].set(xlabel="Terrain Ruggedness Index", ylabel="log GDP (2000)", title="Non African Nations")

ax[1].plot(african_nations["rugged"], african_nations["y_mean"])
ax[1].fill_between(african_nations["rugged"], african_nations["y_perc_5"], african_nations["y_perc_95"], alpha=0.5)
ax[1].plot(african_nations["rugged"], african_nations["true_gdp"], "o")
ax[1].set(xlabel="Terrain Ruggedness Index", ylabel="log GDP (2000)", title="African Nations");

我们观察到，我们的模型和90%置信区间的结果占了我们在实践中观察到的大多数数据点，但仍有相当多的非非洲国家被我们的近似后验概率认为不太可能。

示例:用全秩指南重新审视贝叶斯回归¶

为了改善我们的结果，我们将尝试使用一个指南，从所有参数的多元正态分布中生成样本。这允许我们通过满秩协方差矩阵来捕捉潜在变量之间的相关性Σ∈R5×5；我们之前的指南忽略了这些相关性。也就是说，我们有

α,βa,βr,βar,σu∼qϕ=(μ,Σ)(α,βa,βr,βar,σu)=Normal((α,βa,βr,βar,σu)|μ,Σ)

σ=constrain(σu)

要手动编写这种形式的指南，我们需要组合所有潜在的变量，这样我们就可以pyro.sample他们一起从一个单一的pyro.distributions.MultivariateNormal分发，选择一个实现constrain()确定…的价值σ为积极，创建并初始化参数μ,Σ适当的形状，并限制变化的参数Σ以在整个优化过程中保持有效的协方差矩阵(即保持正定)。

这将是非常乏味的，因此我们将使用另一个自动向导来为我们处理所有这些簿记工作，pyro . infer . auto guide . automultivariatenormal:

[61]:

mvn_guide = pyro.infer.autoguide.AutoMultivariateNormal(model)

使用pyro.render_model证明了不同于我们的平均场AutoNormal指南，该指南明确地捕捉了我们模型中所有潜在变量之间的相关性。新的_AutoMultivariateNormal_latent可视化图形中的节点对应于上面的等式；对应于模型变量的其他节点简单地索引到该张量值随机变量的单个元素中。

[62]:

pyro.render_model(mvn_guide, model_args=(is_cont_africa, ruggedness, log_gdp), render_params=True)

[62]:

我们的模型以及我们的推理和评估代码的其余部分基本上与以前没有变化:我们使用pyro.optim.Adam和pyro.infer.Trace_ELBO以适应新指南的参数，然后从指南中取样并使用Predictive从后验预测分布中取样。

有一个小的区别值得注意:我们通过将引导样本直接传递给Predictive通过itsposterior_samples关键字参数，而不是像上一节那样传递指南。这避免了不必要的重复计算。

[63]:

%%time
pyro.clear_param_store()
mvn_guide = pyro.infer.autoguide.AutoMultivariateNormal(model)
svi = pyro.infer.SVI(model,
                     mvn_guide,
                     pyro.optim.Adam({"lr": 0.02}),
                     pyro.infer.Trace_ELBO())

losses = []
for step in range(1000 if not smoke_test else 2):
    loss = svi.step(is_cont_africa, ruggedness, log_gdp)
    losses.append(loss)
    if step % 100 == 0:
        logging.info("Elbo loss: {}".format(loss))

plt.figure(figsize=(5, 2))
plt.plot(losses)
plt.xlabel("SVI step")
plt.ylabel("ELBO loss")

with pyro.plate("samples", 800, dim=-1):
    mvn_samples = mvn_guide(is_cont_africa, ruggedness)

mvn_gamma_within_africa = mvn_samples["bR"] + mvn_samples["bAR"]
mvn_gamma_outside_africa = mvn_samples["bR"]

# Interface note: reuse guide samples for prediction by passing them to Predictive
# via the posterior_samples keyword argument instead of passing the guide as above
assert "obs" not in mvn_samples
mvn_predictive = pyro.infer.Predictive(model, posterior_samples=mvn_samples)
mvn_predictive_samples = mvn_predictive(is_cont_africa, ruggedness, log_gdp=None)

mvn_gdp = mvn_predictive_samples["obs"]

Elbo loss: 702.4906432628632
Elbo loss: 548.7575962543488
Elbo loss: 490.9642730951309
Elbo loss: 401.81392109394073
Elbo loss: 333.7779414653778
Elbo loss: 247.01823914051056
Elbo loss: 248.3894298672676
Elbo loss: 247.3512134552002
Elbo loss: 248.2095948457718
Elbo loss: 247.21006780862808

CPU times: user 1min 45s, sys: 21.9 ms, total: 1min 45s
Wall time: 7.03 s

现在让我们比较一下前面计算的后验概率AutoDiagonalNormal指南与AutoMultivariateNormal导游。我们将直观地叠加后验分布的横截面(回归系数对的联合分布)。

注意，多元正态近似比平均场近似更分散，并且能够模拟后验系数之间的相关性。

[64]:

svi_samples = {k: v.detach().cpu().numpy() for k, v in samples.items()}
svi_mvn_samples = {k: v.detach().cpu().numpy() for k, v in mvn_samples.items()}

fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(12, 6))
fig.suptitle("Cross-sections of the Posterior Distribution", fontsize=16)
sns.kdeplot(x=svi_samples["bA"], y=svi_samples["bR"], ax=axs[0], bw_adjust=4 )
sns.kdeplot(x=svi_mvn_samples["bA"], y=svi_mvn_samples["bR"], ax=axs[0], shade=True, bw_adjust=4)
axs[0].set(xlabel="bA", ylabel="bR", xlim=(-2.8, -0.9), ylim=(-0.6, 0.2))

sns.kdeplot(x=svi_samples["bR"], y=svi_samples["bAR"], ax=axs[1],bw_adjust=4 )
sns.kdeplot(x=svi_mvn_samples["bR"], y=svi_mvn_samples["bAR"], ax=axs[1], shade=True, bw_adjust=4)
axs[1].set(xlabel="bR", ylabel="bAR", xlim=(-0.55, 0.2), ylim=(-0.15, 0.85))


for label, color in zip(["SVI (Diagonal Normal)", "SVI (Multivariate Normal)"], sns.color_palette()[:2]):
    plt.plot([], [],
                label=label, color=color)

fig.legend(loc='upper right')

[64]:

<matplotlib.legend.Legend at 0x7f8971b854c0>

通过重复我们对非洲内外国家的坚固性-GDP系数分布的可视化，我们可以看到这一点的含义。这两个系数的后验不确定性现在大致相同，与目测数据所暗示的一致。

[65]:

fig, axs = plt.subplots(nrows=1, ncols=2, figsize=(12, 6))
fig.suptitle("Density of Slope : log(GDP) vs. Terrain Ruggedness");

sns.histplot(gamma_within_africa.detach().cpu().numpy(), ax=axs[0], kde=True, stat="density", label="African nations")
sns.histplot(gamma_outside_africa.detach().cpu().numpy(), ax=axs[0], kde=True, stat="density", color="orange", label="Non-African nations")
axs[0].set(title="Mean field", xlabel="Slope of regression line", xlim=(-0.6, 0.6), ylim=(0, 11))

sns.histplot(mvn_gamma_within_africa.detach().cpu().numpy(), ax=axs[1], kde=True, stat="density", label="African nations")
sns.histplot(mvn_gamma_outside_africa.detach().cpu().numpy(), ax=axs[1], kde=True, stat="density", color="orange", label="Non-African nations")
axs[1].set(title="Full rank", xlabel="Slope of regression line", xlim=(-0.6, 0.6), ylim=(0, 11))

handles, labels = axs[1].get_legend_handles_labels()
fig.legend(handles, labels, loc='upper right');

我们对两种近似下非非洲国家的后验预测分布的90%可信区间进行了可视化，验证了我们对观察数据的覆盖有所改善:

[66]:

mvn_predictions = pd.DataFrame({
    "cont_africa": is_cont_africa,
    "rugged": ruggedness,
    "y_mean": mvn_gdp.mean(dim=0).detach().cpu().numpy(),
    "y_perc_5": mvn_gdp.kthvalue(int(len(mvn_gdp) * 0.05), dim=0)[0].detach().cpu().numpy(),
    "y_perc_95": mvn_gdp.kthvalue(int(len(mvn_gdp) * 0.95), dim=0)[0].detach().cpu().numpy(),
    "true_gdp": log_gdp,
})
mvn_non_african_nations = mvn_predictions[mvn_predictions["cont_africa"] == 0].sort_values(by=["rugged"])

fig, ax = plt.subplots(nrows=1, ncols=2, figsize=(12, 6), sharey=True)
fig.suptitle("Posterior predictive distribution with 90% CI", fontsize=16)

ax[0].plot(non_african_nations["rugged"], non_african_nations["y_mean"])
ax[0].fill_between(non_african_nations["rugged"], non_african_nations["y_perc_5"], non_african_nations["y_perc_95"], alpha=0.5)
ax[0].plot(non_african_nations["rugged"], non_african_nations["true_gdp"], "o")
ax[0].set(xlabel="Terrain Ruggedness Index", ylabel="log GDP (2000)", title="Non African Nations: Mean-field")

ax[1].plot(mvn_non_african_nations["rugged"], mvn_non_african_nations["y_mean"])
ax[1].fill_between(mvn_non_african_nations["rugged"], mvn_non_african_nations["y_perc_5"], mvn_non_african_nations["y_perc_95"], alpha=0.5)
ax[1].plot(mvn_non_african_nations["rugged"], mvn_non_african_nations["true_gdp"], "o")
ax[1].set(xlabel="Terrain Ruggedness Index", ylabel="log GDP (2000)", title="Non-African Nations: Full rank");

后续步骤¶

如果你做到了这一步，你就可以开始使用Pyro了！跟随首页上安装Pyro的说明检查一下我们其余的例子和教程，尤其是实用烟火和PyTorch教程系列，包括本教程中相同贝叶斯回归分析的一个版本用一个更py torch-本地建模API.

要了解更多关于Pyro中变分推理的数学背景，请查看我们的SVI教程系列，从第一部分。如果你是PyTorch或深度学习的新手，你也可以从阅读官方介绍中受益“用PyTorch进行深度学习。”

大多数达到这一点的用户也会发现我们的Pyro中张量形状指南必读书。Pyro大量使用了“阵列广播”烘焙到PyTorch和其他数组库中，以并行化模型和推理算法，虽然最初可能很难理解这种行为，但应用直觉和经验法则将大大有助于使您的体验更加顺畅，并避免令人讨厌的形状错误。