Pyro 从入门到出门
有很多同学在入坑 Link Prediction 后私信我概率分布相关的问题。后来我仔细想了想，发现搞机器学习的同学大（就）部（是）分（我）在线性代数、矩阵论、微积分上没什么问题，但是概率论、信息论的知识还是有很大空白的。因此，这里给大家推荐一个现成的库 Pyro。
从 Pytorch 中的概率分布说起
Pytorch 在 torch.distributions 中为我们提供了常见的概率分布，使用这些概率分布我们可以轻松的生成符合预期的随机变量：
m = Bernoulli(torch.tensor([0.3]))
m.sample()
# tensor([ 0.])
# 有 30% 可能性出现 1，70% 可能性出现 0一个概率分布可以由概率分布函数和参数确定，常用的概率分布有以下几种：
正态（高斯）分布
m = Normal(torch.tensor([0.0]), torch.tensor([1.0]))
m.sample()
# tensor([ 0.1046])指数分布
m = Exponential(torch.tensor([1.0]))
m.sample()
# tensor([ 0.1046])二项（伯努利）分布
m = Bernoulli(torch.tensor([0.3]))
m.sample()
# tensor([ 0.])泊松分布
m = Poisson(torch.tensor([4]))
m.sample()
# tensor([ 3.])有些同学可能只会在初始化参数矩阵的时候用到概率分布。随着神经网络的发展，概率模型被吸收和同化，我们越来越多地在神经网络中使用概率分布建模：
它赋予了模型表达随机事件的能力，也在一定程度上带来了优化的困难。
概率计算图的反向传递
在神经网络中，我们想要让模型调整某个表达式的参数，最小化某个目标函数。Hinton 已经告诉我们，这可以通过反向传播（计算图）完成。在随机模型中，我们还想要让模型调整某个概率分布的参数，最小化某个目标函数。
这时，计算图中就包含随机（抽样）节点了，我们可以对反向传播（计算图）进行改造，主要有两种思路：
Score Function
有些抽样无法写出明确的表达式，比如在 0\sim k 中随机选择一个整数。
对于这种情况，我们可以使用 Score Function 方法进行反向传递：在正向传递时，采集一些样本 x ，使用这些样本计算 f(x) （ f 为采样后的操作）作为抽样的得分，并所有样本的平均分(以概率的对数 log(p(x|\theta)) 为权)，以此作为分布的得分；在反向传递时，以得分作为梯度。
在 Pytorch 中，我们只需要调用 sample 就可以了：
probs = policy_network(state)
# Note that this is equivalent to what used to be called multinomial
m = Categorical(probs)
action = m.sample()
next_state, reward = env.step(action)
loss = -m.log_prob(action) * reward
loss.backward()这样做需要的抽样规模较大(batch 较大)，而且在高维空间优化比较困难(方差大)。
Pathwise Dervative
有些抽样可以分解为独立噪声和仿射变换，写出明确的表达式，比如在正态分布 x\sim N(\mu,\sigma^2) 可以写成 x=\mu + \sigma\varepsilon,\varepsilon \sim N(0,1) 。
对于这种情况，我们可以使用 Pathwise Dervative 进行反向传递，使用 Reparameter Trick 将独立噪声分离出来：在正向传播时，生成一个独立噪声 \varepsilon ，使用 \mu 和 \sigma 对 \varepsilon 进行仿射变换，模拟抽样过程；在反向传播时，使用 f(x) 的梯度信息对 \mu 和 \sigma 进行移动，从而达到调整概率分布的目的。
在 Pytorch 中，我们只需要调用 rsample 就可以了：
params = policy_network(state)
m = Normal(*params)
action = m.rsample()
next_state, reward = env.step(action)
loss = -reward
loss.backward()对于某些难以表示的概率分布，这种方法可能效果较差。
对上述两种方法的推导及证明感兴趣的同学可以看 这里 。就实现而言，Pytorch 已经帮助我们隐藏了这些细节。只需要记住，在增强学习中要调用 sample，在 VAE 中要调用 rsample 就可以了。
其他概率编程语言（PPL）
以前，贝叶斯神经网络在神经网络的浪潮中不算瞩目，但是近年来，自动编码器(VAE)和增强学习越来越火，概率神经网络在自然语言处理(NLP)、自动驾驶等领域攻城略地，直接使用 Pytorch 指定复杂的概率模型（除了上面讲到的，还有其他更麻烦的实现）似乎有点捉襟见肘。
在这种情况下，概率编程语言（PPL）的研究越来越火。科研和工程领域都希望 PPL 能够帮助我们直接指定概率事件，不用再关心其实现和优化细节。其中，最出名的两个库是：
选哪个好呢，这就又回到了 TF 男孩和 PT 男孩之间的真理标准问题大讨论。不过说正经的，个人感觉 Pyro 的实现比较丰富，因此推荐各位同学使用 Pyro 上手，Edward 的使用大同小异。（当然，要是只会 TF 的话也可以先学 Edward）
Pyro 封装了什么
作为 GPU 上的 numpy，Pytorch 最擅长的是 Tensor 的管理、各种矩阵运算和反向传播。但是它在推理算法上的实现比较有限。Pyro 利用 Pytorch 在 GPU 上的反向传播，定义了随机计算图的更新方法。
在使用 Pyro 时，我们 不需要手动区分 sample 和 resample 。
# pytorch
normal = torch.distributions.Normal(loc, scale)
x = normal.rsample()
# pyro
x = pyro.sample("my_sample", pyro.distributions.Normal(loc, scale))在贝叶斯推断中，先验分布反应了实验前对总体参数分布的认识，在获得样本后人们会对这个认识进行调整形成后验分布。只要找到后验分布，就相当于找到了样本（输入）和试验结果（输出）之间的规律，能够进行预测等工作了。
由于后验分布的影响因素较为复杂，往往无法直接求得。 我们一般根据先验知识和试验数据设计一个后验概率模型（model），再使用一个较为简单的分布（guide）去逼近它 。
这里我们举一个抛硬币的例子，说明 Pyro 的实现。在抛硬币实验中，我们抛出 10 次硬币，正面向上的次数为 6（试验数据），反面向上的次数为 4：
# 数据集
data = []
for _ in range(6):
data.append(torch.tensor(1.0))
for _ in range(4):
data.append(torch.tensor(0.0))我们都知道，抛出一枚硬币，正面向上的概率为 1/2（先验知识）。但是这次试验的试验结果让我们怀疑这枚硬币密度不均匀，需要调整正面向上概率的认知。综合先验知识和试验数据，我们可以这样 设计 model ：
def model(data):
# 先验知识，正面向上的概率为 1/2
alpha0 = torch.tensor(10.0)
beta0 = torch.tensor(10.0)
f = pyro.sample("latent_fairness", dist.Beta(alpha0, beta0))
# 试验数据，抽样的结果为实际观测值 6/10
for i in range(len(data)):
pyro.sample("obs_{}".format(i), dist.Bernoulli(f), obs=data[i])
# Pyro 希望用户使用 obs 声明观测值这里的 model 比较简单，我们甚至可以手动计算出真实的后验分布模型。但是，在实际应用中，它可能是一个复杂的非线性映射（比如神经网络）。在这种情况下我们往往无法直接求得后验概率模型，需要 定义一个简单的分布（guide） 去逼近后验概率模型：
# Pyro 将训练参数封装为 pyro.param，在反向传播时 Pyro 会帮我们管理这些参数的更新，参数值可以初始化成任何值。
def guide(data):
alpha_q = pyro.param("alpha_q", torch.tensor(15.0),
constraint=constraints.positive)
beta_q = pyro.param("beta_q", torch.tensor(15.0),
constraint=constraints.positive)
pyro.sample("latent_fairness", dist.Beta(alpha_q, beta_q))如何让定义的 简单分布（guide）逼近模型（model） 呢？我们可以考虑最小化每次采样（pyro.sample）中，guide 和 model 的分布差异：
# Pyro 会在计算时对齐 guide 和 model 中的同名的概率分布
# 在求 loss 时，对应的概率分布的差异会被计算
# model
# name == latent_fairness
f = pyro.sample("latent_fairness", dist.Beta(alpha0, beta0))
# guide
# name == latent_fairness
pyro.sample("latent_fairness", dist.Beta(alpha_q, beta_q))这里提醒一下：model 可以比 guide 复杂得多，但是一定不要忘记 guide 中的具名采样一定要在 model 中出现（可以不是同一种采样方式，但是一定要有对应的名称）。衡量两个采样概率分布的差异我们可以使用 KL 散度，要使 guide 逼近 model 只需要最小化所有采样的 KL 散度之和：
\theta'=\mathop{\arg\min}_{\theta}KL(q_{\theta}(z)||p(z|x))\\
这个公式等价于最大化 ELBO：
\theta'=\mathop{\arg\max}_{\theta}\mathbb{E}_q[logp(x,z)-log(q_{\theta}(z))]\\
我们可以使用 Pyro 提供的 ELBO 损失函数，训练过程与 Pytorch 大同小异：
adam_params = {"lr": 0.0005, "betas": (0.90, 0.999)}
optimizer = Adam(adam_params)
svi = SVI(model, guide, optimizer, loss=Trace_ELBO())
for step in range(n_steps):
svi.step(data)贝叶斯神经网络 的实现与上面的例子类似，只是 model 和 guide 的实现稍微复杂一点，是神经网络的形式。这里看一个线性回归的例子：
# 线性回归
class RegressionModel(nn.Module):
def __init__(self, p):
# p = number of features
super(RegressionModel, self).__init__()
self.linear = nn.Linear(p, 1)
self.factor = nn.Parameter(torch.tensor(1.))
def forward(self, x):
return self.linear(x) + (self.factor * x[:, 0] * x[:, 1]).unsqueeze(1)
# 贝叶斯线性回归
def model(is_cont_africa, ruggedness, log_gdp):
a = pyro.sample("a", dist.Normal(8., 1000.))
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.iarange("data", len(ruggedness)):
pyro.sample("obs", dist.Normal(mean, sigma), obs=log_gdp)
def guide(is_cont_africa, ruggedness, log_gdp):
a_loc = pyro.param('a_loc', torch.tensor(0.))
a_scale = pyro.param('a_scale', torch.tensor(1.),
constraint=constraints.positive)
sigma_loc = pyro.param('sigma_loc', torch.tensor(1.),
constraint=constraints.positive)
weights_loc = pyro.param('weights_loc', torch.randn(3))
weights_scale = pyro.param('weights_scale', 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.Normal(sigma_loc, torch.tensor(0.05)))
mean = a + b_a * is_cont_africa + b_r * ruggedness + b_ar * is_cont_africa * ruggedness对比神经网络和贝叶斯神经网络 ，我们会有这样的感觉：
神经网络是一个无类型的编程语言，只能指定执行的流程（网络结构），而无法指定每个变量（训练参数）的类型（实际含义），完全由训练学出；
而贝叶斯神经网络是一个静态类型的编程语言，除了定义执行的流程（网络结构，即 guide），还能指定变量的类型（后验分布表达式，即 model），训练结果相对可解释。
有很多同学不喜欢这样的静态类型语言，没关系，我们可以像 typescript 一样，杂交~ 自动编码器 就是一个很好的例子：
class VAE(nn.Module):
def __init__(self, z_dim=50, hidden_dim=400, use_cuda=False):
super(VAE, self).__init__()
self.encoder = Encoder(z_dim, hidden_dim)
self.decoder = Decoder(z_dim, hidden_dim)
if use_cuda:
self.cuda()
self.use_cuda = use_cuda
self.z_dim = z_dim
# define the model p(x|z)p(z)
def model(self, x):
pyro.module("decoder", self.decoder)
with pyro.plate("data", x.shape[0]):
z_loc = x.new_zeros(torch.Size((x.shape[0], self.z_dim)))
z_scale = x.new_ones(torch.Size((x.shape[0], self.z_dim)))
z = pyro.sample("latent", dist.Normal(z_loc, z_scale).to_event(1))
loc_img = self.decoder.forward(z)
pyro.sample("obs", dist.Bernoulli(loc_img).to_event(1), obs=x.reshape(-1, 784))
def guide(self, x):
pyro.module("encoder", self.encoder)
with pyro.plate("data", x.shape[0]):
z_loc, z_scale = self.encoder.forward(x)
pyro.sample("latent", dist.Normal(z_loc, z_scale).to_event(1))
def reconstruct_img(self, x):
z_loc, z_scale = self.encoder(x)
z = dist.Normal(z_loc, z_scale).sample()
loc_img = self.decoder(z)
return loc_img我们不在意 encoder 和 decoder 中训练参数的概率分布（类型为 any），只考虑隐变量的概率分布（静态类型），并对这一部分定义 model 和 guide。
vae = VAE()
optimizer = Adam({"lr": 1.0e-3})
svi = SVI(vae.model, vae.guide, optimizer, loss=Trace_ELBO())
def train(svi, train_loader, use_cuda=False):
epoch_loss = 0.
for x, _ in train_loader:
if use_cuda: