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

 Heyang Gong 汉化版本的文档 

causalai.github.io/pyro_zh_tutorial/tensor_shapes.html

这个是pyro的官方中文文档
 

其他相关资料 搜索 

pyro.sample("obs"   

或者Pyro

pyro英语带中文字幕的介绍，

用Ubers Pyro介绍概率编程_哔哩哔哩_bilibili

其他资料 视频 教程 基于Pyro和Pytorch的概率编程入门01——简单线性回归的实现

银色狮子巴壁虎

基于Pyro和Pytorch的概率编程入门01——简单线性回归的实现_哔哩哔哩_bilibili

https://www.bilibili.com/video/BV1VF411k7Lm/?spm_id_from=autoNext&vd_source=7d3841cb81fdca8b20a31fa9ee5ec8ef

 基于Pyro和Pytorch的概率编程入门02——pytorch的结构化数据加载与预处理_哔哩哔哩_bilibili

还有pymc3 视频合集

银色狮子巴壁虎的个人空间-银色狮子巴壁虎个人主页-哔哩哔哩视频

概率模块包
Uber的Pyro基于Pytorch

Google的Edward基于TensorFlow

还有一些独立的像

PyMC3,

Stan,

Pomegranate等等

其他代码链接(python2.7 和pyro比较老的版本)：

https://github.com/Rachnog/Deep-Trading/tree/master/bayesian

为了更深入了解概率编程、贝叶斯模型及其应用，我推荐以下资源给大家：

模式识别和机器学习：

http://www.springer.com/us/book/9780387310732

为黑客设计的贝叶斯方法：

https://www.amazon.com/Bayesian-Methods-Hackers-Probabilistic-Addison-Wesley/dp/0133902838

同时推荐以下python库：

PyMC3：

https://github.com/pymc-devs/pymc3

Edward：

http://edwardlib.org/

Pyro：

http://pyro.ai/
 

安装和安装相关程序包

pip install graphviz

pip install pyro-ppl

第一部分 印度小哥的视频教程截图

一个天气预报的预测，

其中有解释 pyro.sample("obs" 采样是什么意思，就其其实就手机正态分布等分布类型的定义

人工智能深度神经网络中名词解释，

mu（均值） 和 sigma（标准差）概率分布参数的常用术语，“weights“（权重）和“bias“（偏置）

`loc`：分布的均值 - `scale`：分布的标准

详细请看 人工智能深度神经网络中名词解释，mu（均值） 和 sigma（标准差）概率分布参数的常用术语，“weights“（权重）和“bias“（偏置）- `loc`：分布的均值 - `scale`：分布的标准-CSDN博客

注意  [2] [3] [4] 其实是同一个意思

第二部分 官方教程 

提示代码要和版本对应上

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

直接指定概率模型可能很麻烦，并且实现它们可能很容易出错。概率编程语言(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. latent random variables \(z\) \(\Longleftrightarrow\) pyro.sample

2. observed random variables \(x\) \(\Longleftrightarrow\) pyro.sample with the obs keyword argument

3. learnable parameters \(\theta\) \(\Longleftrightarrow\) pyro.param

4. plates \(\Longleftrightarrow\) pyro.plate context managers

5. =============  翻译一下 =================

6. 潜在随机变量z ⟺   pyro.sample

7. 观察随机变量x ⟺ 带有obs关键字参数pyro.sample

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

9. 盘子⟺ pyro.plate上下文管理器

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

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

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

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

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

这段文字描述了在 Pyro 这个概率编程框架中如何定义概率模型，以及在最大似然线性回归这个示例模型中各个组成部分的含义。下面是翻译和解释：

1. **潜在随机变量 \( z \)**：在 Pyro 中，潜在随机变量（latent random variables）是通过 `pyro.sample` 函数来定义的。这些变量通常是模型中的隐变量，其值在生成数据的过程中是随机的，但在模型训练时需要被推断出来。

2. **观察随机变量 \( x \)**：观察随机变量（observed random variables）也是通过 `pyro.sample` 函数定义的，但需要使用 `obs` 关键字参数来指定这些变量的观测值。这意味着这些变量的值是已知的，模型需要根据这些观测值来学习参数。

3. **可学习的参数 \( \theta \)**：在 Pyro 中，可学习的参数（learnable parameters）是通过 `pyro.param` 函数来定义的。这些参数是模型在训练过程中需要学习的量，通常对应于模型的权重和偏置。

4. **盘子（plates）**：在 Pyro 中，盘子（plates）是通过 `pyro.plate` 上下文管理器来定义的。它们用于表示数据的批量处理或重复结构，例如在处理多个观测值时。

5. **示例模型：最大似然线性回归**：在最大似然线性回归的示例中，公式 \( \text{mean} = a + b_a \times \text{is\_cont\_africa} + b_r \times \text{ruggedness} + b_{ar} \times \text{is\_cont\_africa} \times \text{ruggedness} \) 表示线性回归模型的预测值（mean），其中：
   - \( a \)：截距项，对应于偏置（bias）。
   - \( b_a \)：与是否位于非洲大陆的特征（is\_cont\_africa）相关的权重。
   - \( b_r \)：与地形崎岖度（ruggedness）相关的权重。
   - \( b_{ar} \)：与是否位于非洲大陆和地形崎岖度交互项相关的权重。

在最大似然线性回归中，我们通常假设数据是独立同分布的，并且通过最小化预测值和观测值之间的差异（通常是平方误差）来估计模型参数。在 Pyro 中，我们可以通过定义这些参数和变量，然后使用相应的概率分布来构建模型，并使用概率推断方法来估计参数。
 

[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)

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

这个报错

如果您想要追踪模型的执行并打印出参数的分布，您可以这样做：

python复制

# 使用 trace 来追踪模型的执行 trace = pyro.poutine.trace(simple_model).get_trace(*model_args) # 打印出参数的分布 print(trace.nodes)

[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 中如何使用 `pyro.sample` 来定义模型中的随机变量，以及如何通过最大似然估计和贝叶斯方法来处理参数估计问题。下面是对这段文字的理解和解释：

1. **`render_distributions = True`**：当设置为 `True` 时，Pyro 将显示参数的约束条件。例如，标准差 `sigma` 必须是非负的，因此约束条件显示为 `"sigma ∈ GreaterThan(lower_bound=0.0)"`。

2. **最大似然估计与贝叶斯方法**：
   - 最大似然估计（MLE）是一种点估计方法，它通过最大化观测数据的似然函数来估计参数。
   - 贝叶斯方法则提供了参数的不确定性估计，它考虑了参数的先验分布和观测数据，通过贝叶斯定理来更新参数的后验分布。

3. **Pyro 中的概率程序**：
   - Pyro 中的概率程序是基于原始概率分布的样本构建的，使用 `pyro.sample` 来标记。
   - `pyro.sample` 可以代表潜在变量或观察变量，具体取决于是否提供了观测值。

4. **`pyro.sample` 的使用**：
   - 当没有提供观测值 `log_gdp` 时，`pyro.sample` 用于定义潜在变量。
   - 当提供了观测值 `log_gdp` 时，`pyro.sample` 用于定义观察变量。

5. **`pyro.sample` 的行为**：
   - 当没有观测值时，`pyro.sample` 相当于调用底层的 `.sample()` 方法来生成样本。
   - 当有观测值时，`pyro.sample` 的行为是返回观测值本身。

6. **贝叶斯推理**：
   - 在 Pyro 中，当模型中的 `pyro.sample` 语句被赋予观测值时，模型的每个其他示例的累积效果都会发生变化，遵循贝叶斯法则。
   - Pyro 的推理算法会“逆向运行程序”，为所有程序分配数学上一致的值。

7. **`pyro.sample` 的命名**：
   - `pyro.sample` 需要一个名称，这是概率编程语言的一个关键特性，它允许用户和 Pyro 的内部机制利用名称来区分模型、观测和推理算法的不同部分。

8. **高阶原语 `pyro.condition`**：
   - `pyro.condition` 是一个高阶原语，它解决了针对单个 Pyro 模型编写多个查询的问题。

总的来说，这段文字强调了 Pyro 中 `pyro.sample` 的灵活性和重要性，以及如何通过贝叶斯方法来处理参数估计中的不确定性。通过命名和条件语句，Pyro 提供了一种强大的方式来构建和推理概率模型。
 

但是，请注意当观察到任何示例语句时，模型中每个其他示例语句的累积效果都会发生变化遵循贝叶斯法则；这是派若的工作推理算法“逆向运行程序”并为所有程序分配数学上一致的值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 中的 pyro.sample 原语的使用，以及如何通过 pyro.condition 高阶原语来处理条件概率模型。下面是对这段文字的解释：

1. 贝叶斯法则和 pyro.sample：

   • 在 Pyro 模型中，当为 pyro.sample 语句提供观测值时，这些观测值会影响模型中其他所有 pyro.sample 语句的累积效应。这是通过贝叶斯法则实现的，即根据新的观测数据更新我们对模型参数的信念。
   • Pyro 的推理算法能够“逆向运行程序”，为模型中的所有 pyro.sample 语句分配数学上一致的值。

2. pyro.sample 的命名：

   • pyro.sample 需要一个名称参数，这是因为名称允许用户和 Pyro 的内部机制区分模型中的不同部分，如潜在变量、观测变量和推理算法。这是概率编程语言的一个关键特性。

3. pyro.condition 高阶原语：

   • pyro.condition 是一个高阶原语，它接受一个模型和一个可能为空的字典，字典中的键是 pyro.sample 语句的名称，值是对应的观测值。
   • 使用 pyro.condition 可以创建一个条件化模型，该模型在给定观测值的情况下进行推理。

4. 使用 pyro.condition 的例子：

   • 假设有一个简单的线性回归模型 simple_model，它接受 is_cont_africa 和 ruggedness 作为输入，并且有一个观测值 log_gdp。
   • 通过 pyro.condition，我们可以移除 log_gdp 作为参数，并用一个更简单的接口 simpler_model 来代替它，如下所示：

     python

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

   • 这里 conditioned_model 是一个条件化模型，它在给定 log_gdp 的情况下进行推理。

5. conditioned_model 与 functools.partial 的比较：

   • 文中提到 conditioned_model 相当于使用 functools.partial 来固定 simple_model 的 log_gdp 参数：

     python

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

   • 然而，这两者之间有一个关键区别：pyro.condition 不仅固定了参数，还告诉 Pyro 这是一个观测值，而 functools.partial 只是创建了一个部分应用的函数。

总结来说，这段文字强调了在 Pyro 中使用 pyro.sample 和 pyro.condition 来构建和推理条件概率模型的重要性。通过命名和条件化，我们可以更灵活地处理模型中的观测数据，并进行有效的贝叶斯推理。

背景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 中的 `pyro.param` 原语及其用法，下面是对这段文字的解读：

1. **`pyro.param` 函数的作用**：
   - `pyro.param` 是 Pyro 中用来读取和写入参数的函数。这些参数存储在 Pyro 的前端键值参数存储区中。

2. **参数定义**：
   - `pyro.param` 的第一个参数是名称（`name`），用于唯一标识参数。
   - 第二个参数是初始值（`init`），可以是一个 `torch.Tensor` 对象或者是一个无参函数，该函数返回一个张量。如果参数已经被初始化，可以省略 `init`。

3. **参数的存储和检索**：
   - 当第一次使用特定名称调用 `pyro.param` 时，它会存储由 `init` 指定的初始值，并返回这个值。
   - 之后，再次使用该名称调用 `pyro.param` 时，它会直接从参数存储中返回该值，而不需要再次指定 `init`。

4. **参数的初始化形式**：
   - 初始值参数 `init` 可以是一个张量对象，也可以是一个返回张量的函数。后者特别有用，因为它可以避免在模型第一次运行时重复构造初始值。

5. **参数的约束**：
   - `pyro.param` 的第三个可选参数 `constraint` 用于指定参数的约束条件。由于许多概率分布的参数具有受限的域，例如正态分布的 `scale` 参数必须是正数，因此约束是一个重要的特性。
   - Pyro 提供了大量预定义的约束，可以在参数初始化时指定，并在每次更新后重新应用这些约束。

6. **参数更新**：
   - 参数存储区可以通过优化算法更新，或者通过调用 `pyro.clear_param_store()` 来清空。
   - 与 `pyro.sample` 不同，`pyro.param` 在模型中可以使用相同的名称多次调用，每次调用都会返回相同的值。

7. **全局参数存储**：
   - 全局参数存储可以通过调用 `pyro.get_param_store()` 来访问。

8. **`pyro.param` 与 PyTorch 的兼容性**：
   - 虽然本教程使用了 `pyro.param` 进行参数管理，但 Pyro 也兼容 PyTorch 的 `torch.nn.Module` API，通过 `pyro.module` 实现。

9. **例子**：
   - 文中还提供了一个使用 `pyro.param` 的示例，展示了如何在 Pyro 模型中声明参数 `a`，并与一个简单的 Python 字典实现进行比较。

总结来说，`pyro.param` 是 Pyro 中用于参数管理和约束的核心原语，它提供了参数的初始化、存储、检索和约束等关键功能，并且与 PyTorch 框架兼容。
 

pyro.plate介绍

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部分，它也有助于管理数据子采样并且作为推断算法的信号，表明某些变量是独立的。

这段文字介绍了 Pyro 中的 `pyro.plate` 原语及其作用，下面是对这段文字的解读：

1. **`pyro.plate` 的定义**：
   - `pyro.plate` 是 Pyro 中用于编码“平板”（plate）符号的原语，用于简化概率机器学习中独立同分布（i.i.d.）随机变量的模型化。

2. **`pyro.plate` 的使用**：
   - `pyro.plate` 可以看作是 Python 中 `for` 循环的替代品，用于在模型中重复执行相同的操作多次，但以一种更高效的方式。

3. **`pyro.plate` 的参数**：
   - `name`：平板的名称，用于区分不同的平板。
   - `size`：平板的大小，即重复执行操作的次数。
   - `dim`：可选参数，指定在输出张量的哪个维度上展开。
   - `other_kwargs`：其他关键字参数。

4. **`pyro.plate` 与 `for` 循环的比较**：
   - 使用 `pyro.plate` 可以替换模型中的 `for` 循环，例如，将以下代码：
     ```python
     with pyro.plate("data", len(ruggedness)):
         return pyro.sample("obs", dist.Normal(mean, sigma), obs=log_gdp)
     ```
     替换为等价的 `for` 循环：
     ```python
     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
     ```

5. **矢量化操作**：
   - 当重复变量的数量很大时，使用 Python 循环会非常慢。`pyro.plate` 利用 PyTorch 的数组广播功能，通过矢量化操作并行执行迭代，从而提高效率。

6. **等价的矢量化代码**：
   - 例如，将 `mean` 和 `sigma` 张量扩展到与 `ruggedness` 长度相同的维度，然后使用 `pyro.sample`：
     ```python
     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)
     ```

7. **`pyro.plate` 的其他用途**：
   - 除了作为矢量化工具外，`pyro.plate` 还有助于管理数据子采样，并作为推断算法的信号，表明某些变量是独立的。

总结来说，`pyro.plate` 是 Pyro 中用于高效处理 i.i.d. 随机变量的原语，它可以替换模型中的 `for` 循环，通过矢量化操作提高效率，并在概率编程中提供额外的语义信息。
 

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

为了使我们的线性回归模型贝叶斯化，我们需要指定先验分布关于参数α∈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中的推理

背景:变分推理¶

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

虽然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中模型和向导（guide）定义的说明。Pyro是一个用于概率编程的库，它允许用户定义复杂的统计模型，并使用各种方法进行推理。在Pyro中，模型（model）和向导（guide）是两种核心组件，它们共同用于变分推断（variational inference）。

下面是对这段文字的解读：

1. **模型定义**：在Pyro中，模型是通过定义一个函数来实现的，这个函数使用`pyro.sample`来声明随机变量。例如，在`model`函数中，`pyro.sample("z_1", ...)`表示声明了一个名为`z_1`的随机变量，并且这个随机变量的分布需要用户指定。

2. **向导定义**：向导也是通过定义一个函数来实现的，这个函数同样使用`pyro.sample`来声明随机变量。在`guide`函数中，`pyro.sample("z_1", ...)`表示声明了一个与模型中同名的随机变量`z_1`。向导的目的是提供一个简化的或者更易于优化的分布来近似模型中的分布。

3. **名称匹配**：模型和向导中的随机变量名称必须一一对应。这意味着在模型中声明的每个随机变量都必须在向导中有相应的声明。

4. **分布的灵活性**：虽然模型和向导可以使用不同的分布来定义相同的随机变量，但是它们提供了很大的灵活性。用户可以根据自己的需求和优化目标来选择合适的分布。

5. **自动制导**：编写向导可能会很困难和乏味，尤其是对于新用户来说。因此，Pyro提供了自动制导的功能，它可以自动生成向导。这包括使用Pyro内置的自动制导方法，或者使用配方自动制导，后者可以自动从模型中生成公共导轨族。

6. **下一节内容预告**：这段文字的结尾提到了下一节将演示自动制导的两种方法，这意味着文档的后续部分将提供具体的示例和指导，帮助用户理解如何使用自动制导。

总的来说，这段文字是对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)

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

进一步阐述了Pyro中的自动向导（autoguide）的概念和作用，以及如何使用它们来进行变分推断。下面是对这段文字的解读：

自动向导集合：Pyro提供了一个广泛的自动向导集合，这些自动向导可以根据给定的模型自动生成向导程序。这意味着用户不需要手动编写向导函数，Pyro可以自动完成这一过程。

AutoGuide实例：pyro.autoguide.AutoGuide是一个自动向导的基类，它本身是一个函数，接受与模型相同的参数，并返回一个包含pyro.sample调用的字典。这个字典定义了向导中的随机变量及其分布。

AutoNormal类：AutoNormal是Pyro提供的一个自动引导类，它能够在一行代码中自动生成一个向导。这个向导相当于用户手动编写的向导，但是更加简单和方便。

向导与推理算法：虽然向导定义了可能的近似后验分布的搜索空间和初始点，但它本身并没有完全指定推理算法。推理算法需要通过优化问题来解决，目的是将初始分布调整到更接近真实后验分布的位置。

优化问题：在变分推断中，我们需要解决一个优化问题，即最大化证据下界（Evidence Lower Bound, ELBO）。这个优化问题涉及到调整向导中的参数，以便使得近似后验分布尽可能地接近模型的真实后验分布。

背景知识：这段文字提到了估计和优化证据下界作为背景知识。证据下界是变分推断中的一个关键概念，它是模型证据的下界，可以通过优化向导参数来提高。

后续内容预告：最后，这段文字预告了接下来的内容将涉及如何制定和解决优化问题，以及如何通过优化来提高变分推断的效果

背景:估计和优化证据下限(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和其他数组库中，以并行化模型和推理算法，虽然最初可能很难理解这种行为，但

应用直觉和经验法则将大大有助于使您的体验更加顺畅，并避免令人讨厌的形状错误。

=====================================================================

其他参考文章 

=========================使用Pytorch和Pyro实现贝叶斯神经网络(Bayesian Neural Network)_gang_akarui-开放原子开发者工作坊============================================

openatomworkshop.csdn.net/664edcdfb12a9d168eb6f5c2.html

第三部分 bnn mnist数据集上的的实现,使用Pytorch和Pyro实现贝叶斯神经网络(Bayesian Neural Network)

如果是 2019年的mnist上面实现的pyro的代码  请把pyro的版本只能装为  0.3.2版本

pip install pyro-ppl==0.3.2

本文翻译自：https://www.noconote.work/entry/2019/01/01/160100最近概率模型和神经网络相结合的研究变得多了起来，这次使用Uber开源的Pyro来实现一个贝叶斯神经网络。概率编程框架最近出了不少，Uber的Pyro基于Pytorch，Google的Edward基于TensorFlow，还有一些独立的像PyMC3,Stan,Pomegranate等等。..

​编辑 gang_akarui  ·  2019-07-18 20:59:17 发布

本文翻译自：https://www.noconote.work/entry/2019/01/01/160100

最近概率模型和神经网络相结合的研究变得多了起来，这次使用Uber开源的Pyro来实现一个贝叶斯神经网络。概率编程框架最近出了不少，Uber的Pyro基于Pytorch，Google的Edward基于TensorFlow，还有一些独立的像PyMC3,Stan,Pomegranate等等。
这次选择使用Pyro实现的原因是，Pyro基于Numpy实现，加上对Pytorch的支持很友好，可以自动计算梯度，动态计算，这些好处当然都是Pytorch带来的。官网链接：https://pyro.ai/， 介绍到此为止，接下来开始一步一步写代码了。

Step1：导入需要的模块

import numpy as np
import torch
import torch.nn as nn
import torch.nn.functional as F
from torchvision import datasets
from torchvision import transforms

import pyro
from pyro.distributions import Normal
from pyro.distributions import Categorical
from pyro.optim import Adam
from pyro.infer import SVI
from pyro.infer import Trace_ELBO

import matplotlib.pyplot as plt
import seaborn as sns
sns.set()

Step2：准备数据集，这次使用的是MNIST

train_loader = torch.utils.data.DataLoader(
        datasets.MNIST("/tmp/mnist", train=True, download=True,
                       transform=transforms.Compose([transforms.ToTensor(),])),
        batch_size=128, shuffle=True)

test_loader = torch.utils.data.DataLoader(
        datasets.MNIST("/tmp/mnist", train=False, transform=transforms.Compose([transforms.ToTensor(),])

Step3：定义神经网络模型

class BNN(nn.Module):
    def __init__(self, input_size, hidden_size, output_size):
        super(BNN, self).__init__()
        self.fc1 = nn.Linear(input_size, hidden_size)
        self.out = nn.Linear(hidden_size, output_size)
        
    def forward(self, x):
        output = self.fc1(x)
        output = F.relu(output)
        output = self.out(output)
        return output

net = BNN(28*28, 1024, 10)

Step4：基于正态分布，初始化weights和bias

def model(x_data, y_data):
    # define prior destributions
    fc1w_prior = Normal(loc=torch.zeros_like(net.fc1.weight),
                                      scale=torch.ones_like(net.fc1.weight))
    fc1b_prior = Normal(loc=torch.zeros_like(net.fc1.bias),
                                      scale=torch.ones_like(net.fc1.bias))
    outw_prior = Normal(loc=torch.zeros_like(net.out.weight),
                                       scale=torch.ones_like(net.out.weight))
    outb_prior = Normal(loc=torch.zeros_like(net.out.bias),
                                       scale=torch.ones_like(net.out.bias))
    
    priors = {
        'fc1.weight': fc1w_prior,
        'fc1.bias': fc1b_prior, 
        'out.weight': outw_prior,
        'out.bias': outb_prior}
    
    lifted_module = pyro.random_module("module", net, priors)
    lifted_reg_model = lifted_module()
    
    lhat = F.log_softmax(lifted_reg_model(x_data))
    pyro.sample("obs", Categorical(logits=lhat), obs=y_data)

Step5：为了近似后验概率分布，先定义一个函数

def guide(x_data, y_data):
    fc1w_mu = torch.randn_like(net.fc1.weight)
    fc1w_sigma = torch.randn_like(net.fc1.weight)
    fc1w_mu_param = pyro.param("fc1w_mu", fc1w_mu)
    fc1w_sigma_param = F.softplus(pyro.param("fc1w_sigma", fc1w_sigma))
    fc1w_prior = Normal(loc=fc1w_mu_param, scale=fc1w_sigma_param)

    fc1b_mu = torch.randn_like(net.fc1.bias)
    fc1b_sigma = torch.randn_like(net.fc1.bias)
    fc1b_mu_param = pyro.param("fc1b_mu", fc1b_mu)
    fc1b_sigma_param = F.softplus(pyro.param("fc1b_sigma", fc1b_sigma))
    fc1b_prior = Normal(loc=fc1b_mu_param, scale=fc1b_sigma_param)

    outw_mu = torch.randn_like(net.out.weight)
    outw_sigma = torch.randn_like(net.out.weight)
    outw_mu_param = pyro.param("outw_mu", outw_mu)
    outw_sigma_param = F.softplus(pyro.param("outw_sigma", outw_sigma))
    outw_prior = Normal(loc=outw_mu_param, scale=outw_sigma_param).independent(1)

    outb_mu = torch.randn_like(net.out.bias)
    outb_sigma = torch.randn_like(net.out.bias)
    outb_mu_param = pyro.param("outb_mu", outb_mu)
    outb_sigma_param = F.softplus(pyro.param("outb_sigma", outb_sigma))
    outb_prior = Normal(loc=outb_mu_param, scale=outb_sigma_param)
    
    priors = {
        'fc1.weight': fc1w_prior,
        'fc1.bias': fc1b_prior,
        'out.weight': outw_prior,
        'out.bias': outb_prior}
    
    lifted_module = pyro.random_module("module", net, priors)
    
    return lifted_module()

Step6：定义Optimizer
*这里多说两句，优化器还是喜闻乐见的Adam；这个SVI函数是干嘛的呢，它使用变分推理（Variational inference）来逼近后验概率分布，说到变分推理自然离不开ELBO，这里的Trace_ELBO函数是继承的父类ELBO，还有其他种类的ELBO，比如JitTrace_ELBO，TraceGraph_ELBO等等，具体选哪个，去看文档就好了。顺便推荐一个变分推理的教程，徐亦达老师的视频教程非常棒。

optim = Adam({"lr": 0.01})
svi = SVI(model, guide, optim, loss=Trace_ELBO())

Step7：训练Bayesian Neural Network

n_iterations = 5
loss = 0

for j in range(n_iterations):
    loss = 0
    for batch_id, data in enumerate(train_loader):
        loss += svi.step(data[0].view(-1,28*28), data[1])

    normalizer_train = len(train_loader.dataset)
    total_epoch_loss_train = loss / normalizer_train
    
    print("Epoch ", j, " Loss ", total_epoch_loss_train)

Step8：使用训练好的模型，预测结果
*Bayesian Neural Network要分别针对每个类别进行采样，然后计算预测值中概率最大的类别，作为结果。这里有10个数字，所以我们设置采样的次数为10

n_samples = 10

def predict(x):
    sampled_models = [guide(None, None) for _ in range(n_samples)]
    yhats = [model(x).data for model in sampled_models]
    mean = torch.mean(torch.stack(yhats), 0)
    return np.argmax(mean.numpy(), axis=1)

print('Prediction when network is forced to predict')
correct = 0
total = 0
for j, data in enumerate(test_loader):
    images, labels = data
    predicted = predict(images.view(-1,28*28))
    total += labels.size(0)
    correct += (predicted == labels).sum().item()
print("accuracy: %d %%" % (100 * correct / total))

*当然也可以把计算出的概率归一化一下，这样使得结果更具有解释性

def predict_prob(x):
    sampled_models = [guide(None, None) for _ in range(n_samples)]
    yhats = [model(x).data for model in sampled_models]
    mean = torch.mean(torch.stack(yhats), 0)
    return mean
#正则化
def normalize(x):
    return (x - x.min()) / (x.max() - x.min())
#可视化
def plot(x, yhats):
    fig, (axL, axR) = plt.subplots(ncols=2, figsize=(8, 4))
    axL.bar(x=[i for i in range(10)], height= F.softmax(torch.Tensor(normalize(yhats.numpy()))[0]))
    axL.set_xticks([i for i in range(10)], [i for i in range(10)])
    axR.imshow(x.numpy()[0])
    plt.show()

效果如下：

x, y = test_loader.dataset[0]
yhats = predict_prob(x.view(-1, 28*28))
print("ground truth: ", y.item())
print("predicted: ", yhats.numpy())
plot(x, yhats)

ground truth:  7
predicted:  [[  39.87391   -86.72453   112.25181   102.81404   -89.95473    52.947186
  -167.88455   402.66757   -29.799454   95.75331 ]]

=========================================================================相关文章2

第四部分 pyro金融方面的应用， 实现贝叶斯神经网络（代码比较老和现在Pyro版本不兼容） 

======================================================

大数据文摘作品

编译：修竹、笪洁琼、夏雅薇

作者用了一种新奇的方法来训练神经网络。更新权重的分布而不是顺序更新静态权重，得到了更有趣和可靠的结果。贝叶斯方法给了我们一个机会，使得我们可以不手动添加正则项的情况下对神经网络进行正则化，理解模型的不确定性，并尽可能使用更少的数据得到更好的结果。

Hi！又见面啦。去年我推出了几篇基于神经网络的金融预测教程，我认为有些结果还是蛮有趣的，值得应用在实际交易中。

如果你读过那些教程，你一定会注意到，当你试图在“随机”数据上用一些机器学习模型并且希望找到隐藏模式时，你其实正逐渐对训练集进行过拟合。

我们使用不同的正则化方法和补充数据来解决这个问题，但是这非常耗时间并且有点盲目搜索了。

今天我想介绍一种稍微不同的方法来用于相同的算法。从概率角度讲，我们可以从数据本身学习正则化方法，在我们预测中估计准确性，使用更少的数据来训练并且在模型中加入概率依赖。

我不会深入到贝叶斯模型或变分推理的技术或者数学细节上，我将给出一些概述，同时也会更加关注如何应用。像往常一样，你可以在下面的链接内查看代码。

代码链接：

https://github.com/Rachnog/Deep-Trading/tree/master/bayesian

为了更深入了解概率编程、贝叶斯模型及其应用，我推荐以下资源给大家：

• 模式识别和机器学习：

  http://www.springer.com/us/book/9780387310732

• 为黑客设计的贝叶斯方法：

  https://www.amazon.com/Bayesian-Methods-Hackers-Probabilistic-Addison-Wesley/dp/0133902838

同时推荐以下python库：

• PyMC3：

  https://github.com/pymc-devs/pymc3

• Edward：

  http://edwardlib.org/

• Pyro：

  http://pyro.ai/

概率编程

这个概率性的东西是什么，而且我们为什么要称之为编程呢？首先，我们先回忆一下“正常”的神经网络以及我们能从中获得什么。

我们有参数（权重），这些参数以矩阵表示，输出通常是一些标量值或者向量（例如用于分类时）。比如说，在用SGD训练模型之后，我们有了这些固定矩阵和网络在相同的输入样本上输出相同的向量。完全正确！

但是如果我们认为这些参数和输出都是互相依赖的分布呢？

神经网络中的每个权重都是来自某个分布的样本，输出也一样，每个输入来自整个网络的样本，同时这个网络依赖参数的样本。它给予了我们什么？

我们从最基础的开始讲。

如果我们把网络看成一组彼此依赖的分布，首先定义联合概率分布为p(y, z|x)，输出为y，还有一些依赖输入x 的模型“内部”、隐藏参数z（和普通神经网络一样）。

而我们想要找到一种神经网络分布，我们可以对y ~ p(y|x)采样然后把分布作为输出（该分布的样本期望值通常就是输出，标准差用来评估不确定性，如果分布模型的尾部越大，我们对于输出越没有信心）。

这个设置或多或少已经很清楚了，我们只需要记住，现在所有的参数，不管是模型的输入还是输出，都是分布。我们需要的训练是找到这些分布的参数以便在实际任务中获得更高的准确率。

必须要提到的是，参数分布的形状是我们自己设置的（例如，所有的初始权重都是w ~ Normal(0, 1)，然后我们将学习正确的均值和方差）。

初始分布称之为先验分布，使用过训练数据拟合参数的分布叫做后验分布。后者用于取样和获得输出数据。

模型的拟合效果怎么样呢？一般的框架叫做变分推理。我们不去深入了解细节，在这里我们需要寻找的模型是可以最大化似然函数log p_w(z|x)的, w 是模型的参数（分布参数），z是隐藏变量（隐藏神经元输出，从参数为 w 的分布中取样得到的），x是输入数据样本。这就是我们的模型。

在Pyro库中我们引入了一个实例作为这个模型的指导，指导中包括一些对所有隐藏变量q_ф(z)的分布，其中 ф叫做变分参数。这个分布必须近似于拟合数据最好的模型参数的“真实”分布。

训练的目的是最小化一个指导中关于输入数据和样本[log(p_w(z|x))—log(q_ф(z))] 的期望。我们这里不讨论训练过程的细节，因为这里面包含好几门大学课程，现在我们就做黑盒优化好了。

哦对了，为什么是编程呢？因为我们通常将这种概率模型（比如神经网络）描述为从一个变量到另一个变量的有向图，这样我们就可以直接表示变量的依赖性：

最初这种概率编程语言被用来定义这些模型并对其进行推断。

为什么用概率编程？

你可以将它认为是一种附加的隐藏变量，从数据中学习模型，而不是采用在模型中注入dropout或L1正则化的方法。

考虑到所有权重都是分布，你可以从中进行N次抽样然后得到输出的分布，通过标准差可以估算你的模型对于结果的准确性。

而且还有个不错的赠礼就是我们只需要用更少的数据来训练模型，并且我们可以在变量间灵活的增加不同的依赖关系。

为什么不用概率编程呢？

我对于使用贝叶斯模型没有太多经验，但就我从Pyro和PyMC3学习中可以知道，训练过程耗时很长而且很难定义准确的先验分布。此外，处理分布的多个样本会导致误解和歧义。

数据展示

我获取了以太坊每日价格。这里面包括典型OHLCV（高开低走）元组以及每天关于以太坊推特的数量。

以太坊价格来源：

https://bitinfocharts.com/

我们将使用7天的价格、交易量和推特数量变化的百分比来预测下一天变化的百分比。

价格、推特数量和交易量变化

上图所示是数据的样本——蓝色表示价格变化，黄色表示推特数量变化，绿色表示交易量变化。这些值（0.1-0.2）之间有一些正相关，所以我们希望利用一些数据训练模型。

贝叶斯线性回归

首先我想了解简单线性回归在我们任务中的表现。

Pyro官方教程：

http://pyro.ai/examples/bayesian_regression.html

我们在PyTorch中定义了我们模型（详细解释请看官方教程）：

class RegressionModel(nn.Module):
def __init__(self, p):
super(RegressionModel, self).__init__()
self.linear = nn.Linear(p, 1)
def forward(self, x):
# x * w + b
return self.linear(x)

这只是我们曾经用过的一个简单确定性模型，但是这也是在Pyro中定义概率的方法。一般线性回归模型分布，均为 ~Normal(0, 1)。我们称之为先验，创建了Pyro的随机函数（在我们例子中，是PyTorch的回归模型），将先验概率加到({‘linear.weight’: w_prior, ‘linear.bias’: b_prior})并且基于输入数据x从模型p(y|x)中采样。

这个模型的指导如下所示：

def model(data):
    # Create unit normal priors over the parameters
    mu = Variable(torch.zeros(1, p)).type_as(data)
    sigma = Variable(torch.ones(1, p)).type_as(data)
    bias_mu = Variable(torch.zeros(1)).type_as(data)
    bias_sigma = Variable(torch.ones(1)).type_as(data)
    
    # Define prior distributions for the weights and bias
    w_prior = Normal(mu, sigma)
    b_prior = Normal(bias_mu, bias_sigma)
    priors = {'linear.weight': w_prior, 'linear.bias': b_prior}
    
    # Lift the regression model with the specified priors
    lifted_module = pyro.random_module("module", regression_model, priors)
    lifted_reg_model = lifted_module()
    
    with pyro.iarange("map", N, subsample=data):
        x_data = data[:, :-1]  # Select the feature columns
        y_data = data[:, -1]   # Select the target column
        
        # Run the regressor forward conditioned on inputs
        prediction_mean = lifted_reg_model(x_data).squeeze()
        
        # Register the observation site with a Normal distribution
        pyro.sample("obs", 
                    Normal(prediction_mean, 
                           Variable(torch.ones(data.size(0))).type_as(data)),
                    obs=y_data.squeeze())

前面是比较老的代码  ，新的代版本可以哟

# 假设 N 是数据的总大小，data 是包含特征和目标列的张量
N = data.size(0)

# 使用 pyro.plate 来声明一个批次，大小为 N，表示我们将会迭代整个数据集
with pyro.plate("map", size=N):
    # 遍历整个数据集
    for i in range(N):
        # 选择当前索引的特征列
        x_data = data[i, :-1]
        
        # 选择当前索引的目标列
        y_data = data[i, -1]
        
        # 运行回归模型，假设 lifted_reg_model 是已经定义好的模型
        prediction_mean = lifted_reg_model(x_data).squeeze()
        
        # 注册观测点，使用 Normal 分布来建模预测和观测之间的分布关系
        # 这里假设我们有一个可学习的观测噪声，用 Variable 来表示
        # 并使用 type_as 来确保它的数据类型与数据张量 data 一致
        observation_noise = Variable(torch.ones(1), requires_grad=True).type_as(data)
        pyro.sample("obs", Normal(prediction_mean, observation_noise), obs=y_data)

上述代码段是使用Pyro库定义一个回归模型的例子。代码中使用了PyTorch作为底层计算框架。下面是格式化后的代码，以及对其中关键部分的解释：

```python

```

关键部分解释：

- `Variable`: 这是Pyro中的一个类，用于创建一个可微分的变量，通常用于参数化模型。
- `torch.zeros` 和 `torch.ones`: 这些是PyTorch函数，用于创建全0或全1的张量。
- `.type_as(data)`: 确保新创建的变量与输入数据`data`的数据类型一致。
- `Normal`: 这是Pyro中定义正态分布的类。
- `priors`: 这是一个字典，定义了模型参数的先验分布。在这个例子中，权重`weight`和偏置`bias`的先验分布都是正态分布。
- `pyro.random_module`: 这个函数用于将一个具有先验分布的模型“提升”为一个Pyro模块。
- `pyro.iarange`: 这是一个上下文管理器，用于指定在哪个索引范围内进行迭代。`"map"`表示全数据集，`N`是数据集的大小，`subsample=data`表示使用的数据子集。
- `lifted_reg_model(x_data)`: 这是模型的前向传播，给定输入`x_data`，模型会计算预测的均值。
- `pyro.sample`: 这个函数用于在Pyro中声明观测数据的分布。这里使用了正态分布，均值为模型预测的均值，方差为1（这里可能需要根据实际情况调整）。

请注意，代码中的`p`, `N`, `regression_model`等变量或函数需要在代码的其他部分定义。此外，`Variable`在Pyro 1.8.0版本之后已经被弃用，建议使用PyTorch的`torch.Tensor`。
 

分布的变分分布。就像你看到的，我们为W和b定义了相同形状的分布，但是尽量使他们更接近实际（只要我们能想到的）。在这个例子里，我选择让这个分布形状更窄一些。

（~Normal(0, 0.1)）

我们以这种方式训练模型：

for j in range(3000):
    epoch_loss = 0.0  # Initialize the epoch loss to zero
    perm = torch.randperm(N)  # Generate a random permutation of indices
    # Shuffle data according to the permutation
    data = data[perm]
    
    # Get indices of each batch
    all_batches = get_batch_indices(N, 64)
    
    for ix, batch_start in enumerate(all_batches[:-1]):  # Exclude the last batch index
        batch_end = all_batches[ix + 1]  # Calculate the end index of the current batch
        batch_data = data[batch_start: batch_end]  # Slice the data for the current batch
        
        # Update the epoch loss with the loss from the current batch
        epoch_loss += svi.step(batch_data)

关键部分解释：

- `for j in range(3000):`: 这是一个循环，迭代3000次，通常用于指定训练模型的轮数（epochs）。
- `epoch_loss = 0.0`: 初始化当前轮的损失为0。
- `torch.randperm(N)`: 生成一个从0到N-1的随机排列，用于数据洗牌。
- `data = data[perm]`: 根据生成的随机排列对数据集进行洗牌。
- `get_batch_indices(N, 64)`: 这个函数用于生成用于批量处理的索引。参数`N`是数据集的大小，64是每个批次的大小。
- `enumerate(all_batches[:-1])`: 通过`enumerate`函数遍历批次索引列表，同时获取索引的序号和值。`[:-1]`确保不包括最后一个批次，这可能是为了避免边界问题。
- `batch_start` 和 `batch_end`: 分别代表当前批次的起始和结束索引。
- `batch_data = data[batch_start: batch_end]`: 根据批次的起始和结束索引，从数据集中提取当前批次的数据。
- `svi.step(batch_data)`: `svi`代表随机变分推断（Stochastic Variational Inference）的实例，`step`方法用于在当前批次上执行一个训练步骤，并返回该批次的损失。

请注意，代码中的`N`, `data`, `get_batch_indices`, `svi`等变量或函数需要在代码的其他部分定义。此外，这段代码假设`svi.step`方法接受批次数据并返回损失值。
 

preds = []

for i in range(100):
    # Sample the regression model parameters from the guide
    sampled_reg_model = guide(X_test)
    
    # Make predictions using the sampled model on the test data
    pred = sampled_reg_model(X_test).data.numpy().flatten()
    
    # Append the flattened prediction to the list of predictions
    preds.append(pred)

之后我们想从模型中取样y。重复取样100次然后计算每一次取样预测的均值和标准差（标准差越大，我们对预测准确的信心越低）。

我们应该记得，金融预测中MSE，MAE或者MAPE等经典指标可能会让人很困惑——相对较小的错误率并不意味着你的模型运行良好，所以在样本外的数据上查看性能是非常重要的，这也是我们要做的：

30天的贝叶斯模型预测

如上图所示，结果并不是很好，但是最后一条的预测形状很好，这给了我们一点信心。让我们继续吧！

普通神经网络

在这个非常简单的模型之后，我们想试着尝试些更有趣的东西，比如神经网络。首先让我们先了解一个简单的MLP，只有一个隐藏层，包括含有25个神经元以及线性激活模型。

下面是格式化后的代码，以及对关键部分的解释：

```python
def get_model(input_size):
    # Define the main input layer with the specified input size
    main_input = Input(shape=(input_size,), name='main_input')
    
    # Add a Dense (fully connected) layer with 25 neurons and linear activation
    x = Dense(25, activation='linear')(main_input)
    
    # Define the output layer with a single neuron and linear activation
    output = Dense(1, activation='linear', name='out')(x)
    
    # Create the Keras model with the input and output layers
    final_model = Model(inputs=[main_input], outputs=[output])
    
    # Compile the model with the 'adam' optimizer and 'mean squared error' loss function
    final_model.compile(optimizer='adam', loss='mse')
    
    return final_model

# Instantiate the model with the size of the input features
model = get_model(len(X_train[0]))

# Train the model for 100 epochs
history = model.fit(
    X_train,              # Training features
    Y_train,              # Training targets
    epochs=100,           # Number of epochs to train for
    batch_size=64,       # Size of each batch during training
    verbose=1,           # Verbosity mode (0 = silent, 1 = progress bar)
    validation_data=(X_test, Y_test),  # Data to validate the model on
    callbacks=[reduce_lr, checkpointer],  # List of callbacks to apply during training
    shuffle=True          # Whether to shuffle the data before each epoch
)
```

关键部分解释：

- `Input`: Keras中的层，用于定义模型的输入。
- `Dense`: Keras中的全连接层，`activation='linear'`指定了线性激活函数。
- `Model`: Keras中的模型类，用于创建模型实例。
- `compile`: 用于配置模型的优化器和损失函数。
- `fit`: 用于训练模型的函数，接受训练数据、验证数据、训练轮数等参数。
- `epochs`: 训练模型的轮数。
- `batch_size`: 每次迭代中使用的样本数量。
- `verbose`: 控制训练过程中的输出信息。
- `validation_data`: 用于模型验证的数据。
- `callbacks`: 训练过程中使用的回调函数列表，例如学习率调整或模型保存。
- `shuffle`: 是否在每个epoch开始前对数据进行洗牌。

请注意，代码中的`X_train`, `Y_train`, `X_test`, `Y_test`, `reduce_lr`, `checkpointer`等变量或对象需要在代码的其他部分定义。此外，`reduce_lr`和`checkpointer`应该是Keras的回调函数实例，用于在训练过程中调整学习率或保存模型。
 

得到以下结果：

30天的Keras神经网络预测

这个结果甚至比简单贝叶斯回归还糟糕，而且这个模型不能得到确定性估计，更重要的是，这个模型甚至不能正则化。

贝叶斯神经网络

现在我想在PyTorch中定义一个和我们在Keras中训练的相同的神经网络。

class Net(torch.nn.Module):
    def __init__(self, n_feature, n_hidden):
        super(Net, self).__init__()
        # Define the hidden layer with linear activation
        self.hidden = torch.nn.Linear(n_feature, n_hidden)
        
        # Define the output layer with a single output unit
        self.predict = torch.nn.Linear(n_hidden, 1)
    
    def forward(self, x):
        # Pass the input through the hidden layer
        x = self.hidden(x)
        
        # Pass the output of the hidden layer through the output layer
        x = self.predict(x)
        
        return x

关键部分解释：

• class Net(torch.nn.Module): 定义了一个名为Net的类，它继承自torch.nn.Module，这是所有神经网络模块的基类。
• def __init__(self, n_feature, n_hidden): 类构造函数，接受输入特征的数量n_feature和隐藏层的神经元数量n_hidden。
• super(Net, self).__init__(): 调用基类的构造函数。
• self.hidden: 定义一个线性层作为网络的隐藏层，输入特征数量为n_feature，输出特征数量为n_hidden。
• self.predict: 定义一个线性层作为网络的输出层，输入特征数量为n_hidden，输出特征数量为1。
• def forward(self, x): 定义前向传播函数，它接受输入x。
• x = self.hidden(x): 将输入x传递到隐藏层。
• x = self.predict(x): 将隐藏层的输出传递到输出层。
• return x: 返回网络的最终输出。

请注意，这个网络模型是一个简单的线性神经网络，没有激活函数。在实际应用中，可能需要在隐藏层添加非线性激活函数以提高模型的表达能力。此外，n_feature应该是输入数据的特征数量，n_hidden是隐藏层的神经元数量。

与贝叶斯回归模型相比，我们现在有两组参数（从输入到隐藏层以及从隐藏层到输出），所以我们稍微改变一下模型分布和先验：

priors = {'hidden.weight': w_prior,
'hidden.bias': b_prior,
'predict.weight': w_prior2,
'predict.bias': b_prior2}

以及这个指导：

dists = {'hidden.weight': w_dist,
'hidden.bias': b_dist,
'predict.weight': w_dist2,
'predict.bias': b_dist2}

不要忘记为模型中所有的分布设置不同的名字，因为不能有任何的歧义和重复！可以在源代码中查看更多细节。

源代码：

https://github.com/Rachnog/Deep-Trading/tree/master/bayesian

在拟合模型和采样后，让我们直接看最终结果：

30天的Pyro神经网络预测

这个结果看上去比之前的结果都要好！

考虑下从贝叶斯模型中学到的正则化或者权重的性质，与普通神经网络做比较，我还会看一下权重统计。下面是我在Pryo模型中如何检查参数的：

for name in pyro.get_param_store().get_all_param_names():
print name, pyro.param(name).data.numpy()

在Keras模型中我是这么做的：

import tensorflow as tf
sess = tf.Session()
with sess.as_default():
tf.global_variables_initializer().run()
dense_weights, out_weights = None, None
with sess.as_default():
for layer in model.layers:
if len(layer.weights) > 0:
weights = layer.get_weights()
if 'dense' in layer.name:
dense_weights = layer.weights[0].eval()
if 'out' in layer.name:
out_weights = layer.weights[0].eval()

举个例子，对于Keras模型中，最后一层的权重平均值和标准差分别是-0.0025901748，0.30395034，对于Pyro模型分别是0.0005974418和0.0005974418，数值更小，模型性能更好！

就像许多L2或者dropout这种正则化方法做的那样，让参数尽可能接近0，然后我们可以用变分推理来实现！对于隐藏层权重的情况就更有趣了。

我们把一些权重向量画出来，蓝色代表Keras的权重，橙色代表Pyro的权重：

输入和隐藏层间的一些权重

有趣的是，事实上不仅权重的均值和标准差很小，而且权重变得更加稀疏，所以基本上我们对于第一组权重用到了稀疏表示（类似L1正则），对第二组用到了类似L2正则表示，简直不可思议！不要忘记试一下代码哦！

结论

我们用了一种新奇的方法来训练神经网络。我们更新权重的分布而不是顺序更新静态权重。所以我们得到更有趣和可靠的结果。

我想要强调的是，贝叶斯方法给了我们一个机会，使得我们可以不手动添加正则项的情况下对神经网络进行正则化，理解模型的不确定性，并尽可能使用更少的数据得到更好的结果。欢迎继续关注!:)

原文链接：

https://medium.com/@alexrachnog/financial-forecasting-with-probabilistic-programming-and-pyro-db68ab1a1dba

手把手：基于概率编程Pyro的金融预测，让正则化结果更有趣！_网易订阅

https://zhuanlan.zhihu.com/p/563715475

其他啊文章baidu 搜索 pypro 或者 

变分推断之变分自编码器（VAE）续 （Pyro 实现）

​

丁贵金：变分推断之变分自编码器（VAE）续 （ELBO）6 赞同 · 0 评论文章​编辑

丁贵金：变分推断之变分自编码器（VAE）续 LOSS FUNCTION8 赞同 · 5 评论文章​编辑

丁贵金：变分推断之变分自编码器（VAE）12 赞同 · 0 评论文章​编辑

做了这么多准备，就是为了能够在 Pyro 中试一试。Pyro 自称的深度概率编程，很大程度上是因为能够和 PyTorch 的神经网络结合使用。

Pyro 本身直接提供了 ELBO 的损失函数 （LOSS FUNCTION），可以省去了写损失函数的部分。

上代码，网络结构基本一致，不是全连接就是卷积，这不是我们关注的重点。

# define the PyTorch module that parameterizes the
# diagonal gaussian distribution q(z|x)
class Encoder(nn.Module):
    def __init__(self, z_dim, hidden_dim):
        super().__init__()
        # setup the three linear transformations used
        self.fc1 = nn.Linear(784, hidden_dim)
        self.fc21 = nn.Linear(hidden_dim, z_dim)
        self.fc22 = nn.Linear(hidden_dim, z_dim)
        # setup the non-linearities
        self.softplus = nn.Softplus()

    def forward(self, x):
        # define the forward computation on the image x
        # first shape the mini-batch to have pixels in the rightmost dimension
        x = x.reshape(-1, 784)
        # then compute the hidden units
        hidden = self.softplus(self.fc1(x))
        # then return a mean vector and a (positive) square root covariance
        # each of size batch_size x z_dim
        z_loc = self.fc21(hidden)
        z_scale = torch.exp(self.fc22(hidden))
        return z_loc, z_scale


# define the PyTorch module that parameterizes the
# observation likelihood p(x|z)
class Decoder(nn.Module):
    def __init__(self, z_dim, hidden_dim):
        super().__init__()
        # setup the two linear transformations used
        self.fc1 = nn.Linear(z_dim, hidden_dim)
        self.fc21 = nn.Linear(hidden_dim, 784)
        # setup the non-linearities
        self.softplus = nn.Softplus()

    def forward(self, z):
        # define the forward computation on the latent z
        # first compute the hidden units
        hidden = self.softplus(self.fc1(z))
        # return the parameter for the output Bernoulli
        # each is of size batch_size x 784
        loc_img = torch.sigmoid(self.fc21(hidden))
        return loc_img

主要看 model 和 guide 函数是怎么写的：

# define a PyTorch module for the VAE
class VAE(nn.Module):
    # by default our latent space is 50-dimensional
    # and we use 400 hidden units
    def __init__(self, z_dim=50, hidden_dim=400, use_cuda=False):
        super().__init__()
        # create the encoder and decoder networks
        self.encoder = Encoder(z_dim, hidden_dim)
        self.decoder = Decoder(z_dim, hidden_dim)

        if use_cuda:
            # calling cuda() here will put all the parameters of
            # the encoder and decoder networks into gpu memory
            self.cuda()
        self.use_cuda = use_cuda
        self.z_dim = z_dim

    # define the model p(x|z)p(z)
    def model(self, x):
        # register PyTorch module `decoder` with Pyro
        pyro.module("decoder", self.decoder)
        with pyro.plate("data", x.shape[0]):
            # setup hyperparameters for prior p(z)
            z_loc = torch.zeros(x.shape[0], self.z_dim, dtype=x.dtype, device=x.device)
            z_scale = torch.ones(x.shape[0], self.z_dim, dtype=x.dtype, device=x.device)
            # sample from prior (value will be sampled by guide when computing the ELBO)
            z = pyro.sample("latent", dist.Normal(z_loc, z_scale).to_event(1))
            # decode the latent code z
            loc_img = self.decoder.forward(z)
            # score against actual images (with relaxed Bernoulli values)
            pyro.sample(
                "obs",
                dist.Bernoulli(loc_img, validate_args=False).to_event(1),
                obs=x.reshape(-1, 784),
            )
            # return the loc so we can visualize it later
            return loc_img

    # define the guide (i.e. variational distribution) q(z|x)
    def guide(self, x):
        # register PyTorch module `encoder` with Pyro
        pyro.module("encoder", self.encoder)
        with pyro.plate("data", x.shape[0]):
            # use the encoder to get the parameters used to define q(z|x)
            z_loc, z_scale = self.encoder.forward(x)
            # sample the latent code z
            pyro.sample("latent", dist.Normal(z_loc, z_scale).to_event(1))

    # define a helper function for reconstructing images
    def reconstruct_img(self, x):
        # encode image x
        z_loc, z_scale = self.encoder(x)
        # sample in latent space
        z = dist.Normal(z_loc, z_scale).sample()
        # decode the image (note we don't sample in image space)
        loc_img = self.decoder(z)
        return loc_img

这里直接封装了一个 VAE 的类，除了 model 和 guide 还提供了一个 reconstruct_img，通过图片生成隐变量，再生成图片。

这里使用了 pyro.module 相当于一组参数 params 注册给 Pyro，这些参数就是神经网络的参数。神经网络的参数是全局的，能够作用在每个数据上，无论在 guide 和 model 中都设置了盘子来控制隐变量的上下文，对于数据，每个图像和隐变量都一一对应。

和以往我们讨论的不一样，guide 中使用了观察数据，因为要借助它来采样，观察数据和encoder网络的参数一起提供隐变量的loc和scale（均值和方差），然后通过 guide 采样 ‘latent' 隐变量。

训练（推理）时候使用了 mini batch，所以通过数据维度x.shape[0]抓到 mini batch 有多少个图片一次性传进来。

to_event(1) 是把第1个维度作为描述事件的维度。因为同时还有batch的维度需要分开。还有一点值得注意的是这decoder输出的是概率，为了采样黑白图，直接用 Bernouli 二项分布。

def main(args):
    # clear param store
    pyro.clear_param_store()

    # setup MNIST data loaders
    # train_loader, test_loader
    train_loader, test_loader = setup_data_loaders(
        MNIST, use_cuda=args.cuda, batch_size=256
    )

    # setup the VAE
    vae = VAE(use_cuda=args.cuda)

    # setup the optimizer
    adam_args = {"lr": args.learning_rate}
    optimizer = Adam(adam_args)

    # setup the inference algorithm
    elbo = JitTrace_ELBO() if args.jit else Trace_ELBO()
    svi = SVI(vae.model, vae.guide, optimizer, loss=elbo)

    # setup visdom for visualization
    if args.visdom_flag:
        vis = visdom.Visdom()

    train_elbo = []
    test_elbo = []
    # training loop
    for epoch in range(args.num_epochs):
        # initialize loss accumulator
        epoch_loss = 0.0
        # do a training epoch over each mini-batch x returned
        # by the data loader
        for x, _ in train_loader:
            # if on GPU put mini-batch into CUDA memory
            if args.cuda:
                x = x.cuda()
            # do ELBO gradient and accumulate loss
            epoch_loss += svi.step(x)

        # report training diagnostics
        normalizer_train = len(train_loader.dataset)
        total_epoch_loss_train = epoch_loss / normalizer_train
        train_elbo.append(total_epoch_loss_train)
        print(
            "[epoch %03d]  average training loss: %.4f"
            % (epoch, total_epoch_loss_train)
        )

        if epoch % args.test_frequency == 0:
            # initialize loss accumulator
            test_loss = 0.0
            # compute the loss over the entire test set
            for i, (x, _) in enumerate(test_loader):
                # if on GPU put mini-batch into CUDA memory
                if args.cuda:
                    x = x.cuda()
                # compute ELBO estimate and accumulate loss
                test_loss += svi.evaluate_loss(x)

                # pick three random test images from the first mini-batch and
                # visualize how well we're reconstructing them
                if i == 0:
                    if args.visdom_flag:
                        plot_vae_samples(vae, vis)
                        reco_indices = np.random.randint(0, x.shape[0], 3)
                        for index in reco_indices:
                            test_img = x[index, :]
                            reco_img = vae.reconstruct_img(test_img)
                            vis.image(
                                test_img.reshape(28, 28).detach().cpu().numpy(),
                                opts={"caption": "test image"},
                            )
                            vis.image(
                                reco_img.reshape(28, 28).detach().cpu().numpy(),
                                opts={"caption": "reconstructed image"},
                            )

            # report test diagnostics
            normalizer_test = len(test_loader.dataset)
            total_epoch_loss_test = test_loss / normalizer_test
            test_elbo.append(total_epoch_loss_test)
            print(
                "[epoch %03d]  average test loss: %.4f" % (epoch, total_epoch_loss_test)
            )

        if epoch == args.tsne_iter:
            mnist_test_tsne(vae=vae, test_loader=test_loader)
            plot_llk(np.array(train_elbo), np.array(test_elbo))

    return vae

主函数的过程没有什么特别，JIT 是用来产生TRACE，然后可以给 C++/C 的程序加载模型用的。

从 Pyro 实现的 VAE 来看，它更清晰的描述了隐变量和数据一一对应的关系，以及神经网络的参数是全局都要依赖的关系（所有数据），隐去了显式计算 ELBO 的过程，但是 guide 和 model 要做好准备。guide 采样的 z 依赖了数据 x，这个和我们以前分析的都不太一样，很多guide是不依赖于数据的，直接给出隐变量的采样，所以 VAE 的 z 实际上是 z｜x （给定 x 的 z）。z 的先验选择了标准正态分布。 model 和 guide 最终完成的事情还是一致的，guide 负责采样隐变量，model 提供先验和对观察数据采样和数据生成。另外 deocder 输出的是概率分布，能够直接用来绘图，因为是灰度图像可以这样做。VAE 从优化 �(�′|�) 的入手，而实现对 �(�|�) 和 �(�′|�) 进行优化，最后将参数落实在神经网络的参数上。