TowardsDataScience 博客中文翻译 2020（一百四十九）

原文：TowardsDataScience Blog

协议：CC BY-NC-SA 4.0

贝叶斯推理——直觉和例子

原文：https://towardsdatascience.com/bayesian-inference-intuition-and-example-148fd8fb95d6?source=collection_archive---------0-----------------------

使用 Python 代码

为什么有人要发明贝叶斯推理？

一句话:随着我们收集更多的数据，更新概率 。

贝叶斯推理的核心是将两个不同的分布(似然和先验)组合成一个“更聪明”的分布(后验)。在经典的最大似然估计(MLE)没有考虑先验的意义上，后验是**“更聪明”。一旦我们计算出后验概率，我们就用它来寻找“最佳”参数，而“最佳”就是在给定数据的情况下，最大化后验概率****。这个过程叫做最大后验概率(MAP) 。MAP 中使用的优化和典型机器学习中使用的优化是一样的，比如梯度下降或者牛顿法等。**

当我研究著名的贝叶斯法则时，通过分析理解这个等式是相当容易的。但是你如何用数据实现它呢？

统计世界中最著名的方程式之一

具体来说，我们会有大量的数据点 X，我们如何将概率 wrt X 乘以概率 wrt θ？****

显然，贝叶斯推理的艺术在于你如何实现它。

示例:

每天大约有 2000 名读者访问我的中型博客。有人看完文章鼓掌有人不鼓掌。我想预测一下，当我将来写一篇新的博文时，会有百分之多少的人参与进来并鼓掌**。**

这类问题适用面广。尝试将这一点应用到你自己的建模工作中——广告的点击率、在你的网站上实际购买的客户的转化率、同意与你约会的人的百分比、患乳腺癌的女性的生存机会等。

生成数据:

让我们生成数据 X 。现实生活中，你对 X 没有任何控制这是你要观察的。

import numpy as np
np.set_printoptions(threshold=100)# Generating 2,000 readers' reponse. 
# Assuming the claps follow a Bernoulli process - a sequence of binary (success/failure) random variables.
# 1 means clap. 0 means no clap.**# We pick the success rate of 30%.**
**clap_prob = 0.3**# IID (independent and identically distributed) assumption
**clap_data = np.random.binomial(n=1, p=clap_prob, size=2000)**

拍手数据是二进制的。1 表示拍手，0 表示没有拍手。

数据 X 看起来像:

In [1]: clap_data
Out[1]: array([0, 0, 0, ..., 0, 1, 0])In [2]: len(clap_data)
Out[2]: 2000

贝叶斯推理有三个步骤。

第一步。【先验】选择一个 PDF 来建模你的参数 θ ，又名先验分布 P(θ) 。这是你在** 看到数据 X 之前关于参数 的最佳猜测。**

第二步。【可能性】为 P(X|θ) 选择一个 PDF** 。基本上你是在给定参数 θ 的情况下模拟数据 X 的样子。**

第三步。【后验】计算后验**分布 P(θ|X) 取 P(θ|X)最高的 θ。**

后验成为新的先验。获得更多数据后，重复第 3 步。

第一步。先验 P(θ)

第一步是选择 PDF** 对参数 θ 建模。**

参数θ代表什么？

鼓掌声👏概率。

那么，我们应该用什么样的概率分布来建模一个概率呢？

要表示一个概率，需要满足几个条件。首先，域的范围应该是从 0 到 1。第二，应该是连续分布。

那么我能想到的有两个众所周知的概率分布:

贝塔和狄利克雷。

狄利克雷用于多变量，贝塔用于单变量。我们只有一件事可以预测，那就是一个概率，所以我们用贝塔分布。

(一个有趣的旁注:很容易在(0，1)上创建任意分布。取任何一个在 0 到 1 之间不爆炸的函数，保持正的。然后，简单的从 0 到 1 积分，用那个结果除函数。)

要使用 Beta 分布，我们需要确定两个参数 α & β 。你可以把 α 想成有多少人鼓掌(成功数)β想成有多少人不鼓掌(失败数)。这些参数——αα&β的大小——将决定分布的形状。

假设您有昨天的数据，并观察到 2000 名访客中有 400 人鼓掌。

根据 beta 分布你会怎么写呢？

import scipy.stats as stats
import matplotlib.pyplot as plt**a = 400
b = 2000 - a**# domain θ
theta_range = np.linspace(0, 1, 1000)# prior distribution P(θ)
prior = **stats.beta.pdf(x = theta_range, a=a, b=b)**

让我们画出所有 θ 值的先验分布。

# Plotting the prior distribution
plt.rcParams['figure.figsize'] = [20, 7]
fig, ax = plt.subplots()
**plt.plot(theta_range, prior,** linewidth=3, color='palegreen')# Add a title
plt.title('[Prior] PDF of "Probability of Claps"', fontsize=20)# Add X and y Label
plt.xlabel('θ', fontsize=16)
plt.ylabel('Density', fontsize=16)# Add a grid
plt.grid(alpha=.4, linestyle='--')# Show the plot
plt.show()

不熟悉 Y 轴上的术语“密度”？→ 阅读“PDF 不是概率”

正如预期的那样，它的峰值达到了 20% (400 次鼓掌/ 2000 个读者)。两千个数据点似乎产生了一个强有力的先验。如果我们使用更少的数据点**，比如说 100 个阅读器，曲线将会不那么尖锐。用α= 20&β= 80 试试。**

对于那些想知道概率密度怎么会大于 1 的人。👉概率密度不是概率。

第二步。可能性 P(X|θ)

选择一个概率模型为 P(X|θ) ，给定特定参数 θ 看到数据的概率 X 。可能性也称为抽样分布。对我来说，“抽样分布”这个术语比“可能性”更直观。

要选择使用哪种概率分布来模拟抽样分布，我们首先需要问:

我们的数据 X 是什么样子的？

X 是二进制数组[0,1,0,1,...,0,0,0,1]。

我们还有访客总数(n)** 和我们想要**鼓掌的概率§。****

好的，n&p……他们对你尖叫什么？

n&p的二项分布。

**# The sampling dist P(X|θ) with a prior θlikelihood = **stats.binom.pmf(k = np.sum(clap_data), n = len(clap_data), p = a/(a+b))****

我们之前的假设 θ 可能性大吗？

**In [63]: likelihoodOut[63]: 4.902953768848812e-30**

没有。

我们来看一下 P(X|θ) 对于所有可能的 θ 的曲线图。

**# Likelihood P(X|θ) **for all θ's**
likelihood = stats.binom.pmf(k = np.sum(clap_data), n = len(clap_data), **p = theta_range**)# Create the plot
fig, ax = plt.subplots()
**plt.plot(theta_range, likelihood,** linewidth=3, color='yellowgreen')# Add a title
plt.title('[Likelihood] Probability of Claps' , fontsize=20)# Add X and y Label
plt.xlabel(’θ’, fontsize=16)
plt.ylabel(’Probability’, fontsize=16)# Add a grid
plt.grid(alpha=.4, linestyle='--')# Show the plot
plt.show()**

给定数据的可能性 P(X |θ)

第三步。后。P(θ|X)

最后，让我们来回答我们在开始时提出的问题:

具体来说，我们会有大量的数据点 X，我们如何将概率 wrt X 乘以概率 wrt θ？

即使有数千个数据点，我们也可以通过将数据插入您选择的模型(在本例中，二项式分布)中，将它们转换为单个标量——可能性 P(X|θ) — 。)

然后，我们对一个具体的 θ 计算P(θ)&P(X |θ)，并将它们相乘。如果你对每一个可能的 θ 都这样做，你可以在不同的 θ 的中挑选最高的 P(θ) * P(X|θ) 。

你最初对参数的猜测是 P(θ)。现在你正在将简单的 P(θ)升级为更具信息性的东西——P(θ| X)——随着更多数据的可用。
P(θ|X) 仍然是 θ 的概率，就像 P(θ) 是一样。不过， P(θ|X) 是 P(θ) 更聪明的版本。

代码胜过千言万语:

# (cont.)theta_range_e = theta_range + 0.001 prior = **stats.beta.cdf(x = theta_range_e, a=a, b=b) - stats.beta.cdf(x = theta_range, a=a, b=b)** # prior = stats.beta.pdf(x = theta_range, a=a, b=b)likelihood = **stats.binom.pmf(k = np.sum(clap_data), n = len(clap_data), p = theta_range)** posterior = **likelihood * prior** # element-wise multiplication
normalized_posterior = posterior / np.sum(posterior)

In [74]: np.argmax(prior)
Out[74]: 199In [75]: np.argmax(likelihood)
Out[75]: 306In [76]: np.argmax(posterior)
Out[76]: 253

# Plotting all three together
fig, axes = plt.subplots(3, 1, sharex=True, figsize=(20,7))
plt.xlabel('θ', fontsize=24)**axes[0].plot(theta_range, prior,** label="**Prior**", linewidth=3, color='palegreen')
axes[0].set_title("Prior", fontsize=16)**axes[1].plot(theta_range, likelihood**, label="**Likelihood**", linewidth=3, color='yellowgreen')
axes[1].set_title("Sampling (Likelihood)", fontsize=16)**axes[2].plot(theta_range, posterior,** label='**Posterior**', linewidth=3, color='olivedrab')
axes[2].set_title("Posterior", fontsize=16)
plt.show()

当你看后验图(第三张)时，注意这是可能性向先验转移的地方。 先前的鼓掌概率为 20%。数据的鼓掌概率被给定为 30%。现在，后验概率的峰值在 0.25%左右。

此外，请注意，在先验/可能性中，钟形曲线的宽度在后验中缩小了。因为我们通过采样整合了更多的信息，所以可能的参数范围现在变窄了。

你收集的数据越多，后验概率图看起来就越像概率图，越不像先验概率图。换句话说，当你得到更多的数据时，原来的先验分布就不那么重要了。

最后，我们选择θ，其给出通过数值优化计算的最高后验**，例如梯度下降或牛顿法。这整个迭代过程被称为最大后验估计(MAP) 。**

脚注:我们通过减去两个stats.beta.**cdf**而不是使用stats.beta.**pdf**来计算先验，因为可能性stats.binom.**pmf** 是一个概率，而stats.beta.**pdf**返回一个密度**。即使我们用密度来计算后验概率，也不会改变优化结果。然而，如果你想要单位匹配，将密度转换成概率是必要的。(还不清楚？阅读“PDF 不是概率，而 PMF 是概率”。)**

需要注意一些事情:

1. MAP 不仅可以通过数值优化来计算，还可以通过(仅适用于某些分布)、修改的期望最大化算法或蒙特卡罗方法来计算。****
2. 使用封闭形式计算后验概率非常方便，因为你不必进行昂贵的计算。(看贝叶斯公式分母中每一个可能的 θ 的积分。非常贵。)
   在我们上面的例子中，贝塔分布是二项似然的共轭先验。这意味着，在建模阶段，我们已经知道后验也将是贝塔分布。那么在进行更多的实验后，我们只需将成功和失败的次数分别加到现有的参数 α 和 β 上，就可以计算出后验。
   我们什么时候可以将这种轻松的封闭形式用于后部？
   当你的先验有一个封闭形式(当然)并且是一个共轭先验——当先验乘以似然性再除以归一化常数，仍然是同一个分布。
3. 归一化(对分母中每个可能的 θ 的积分)对于寻找最大后验概率是否必要？****

没有。不用归一化还是能找到最大值。但是，如果您想比较不同模型的后验概率，或者计算点估计，您需要将其归一化。

4.MAP 估计可以被看作是通过先验引入更多信息的正则化 ML。

5.你熟悉最大似然估计(MLE)，但不太熟悉最大后验概率(MAP)吗？
MLE 只是具有统一先验的映射的一个特例。

6.我写了先验是你最好的猜测关于参数 在 看到数据之前，然而在实际操作中，一旦我们计算了后验，后验就变成了新的先验，直到新的一批数据进来。这样，我们可以迭代地更新我们的先验和后验。马尔可夫链蒙特卡罗抽样方法就是基于这种思想。****

模拟 NBA 选秀顺序的贝叶斯推理

原文：https://towardsdatascience.com/bayesian-inference-of-mock-nba-draft-order-f4fa4a26be17?source=collection_archive---------52-----------------------

如何使用优步概率编程库 Pyro 对 NBA 选秀顺序进行建模

我们许多密切关注 NBA 选秀的人都渴望知道选秀的准确顺序，例如每个球员将被带到哪里(“哪里”在这里有双重含义，球队和选秀位置)。为此，我们通常会从十几份甚至更多的“重要”模拟草稿中吸取信息。分析师和媒体试图以不同的方式将这些信息整合成我喜欢称之为“元模拟”的东西。一般来说，这只涉及非常简单的操作，如平均排名或查看最小和最大吃水位置。Chris Feller 前几天分享了一个非常酷的生存分析方法。我还没有深入研究它，但是和我在这里做的事情相比会很有趣。

我想做的是为模拟草稿创建一个生成元模型，这与使用生存分析略有不同，但可以产生类似的见解，我将在后面展示。理论上，这使得人们能够从等级排序中推断出一些“价值”的度量(我没有发现模拟选秀给玩家分配数字评级)，然后从后验分布中生成随机模拟选秀顺序。换句话说，就像 Nate Silver 运行模拟来预测谁将赢得选举一样，我们可以为模拟草案创建一个模拟，并在其中做各种各样的事情。让我告诉你怎么做！

我们首先需要的是对等级数据有意义的概率分布。向 Plackett-Luce 分布问好。基本思想是，如果你有一些集合 S 的项目，每个项目都有一些“价值” w_i 你可以在一些向量中的这些项目的任意排列的排名上生成一个概率分布。我不会向您展示这些方程(参见参考资料)，但我可以向您展示一个非常简单的模拟来快速理解这个想法。首先我们需要导入一些包并定义 PL 发行版，我发现有人已经在 github 这里上实现了。这里的工作都是使用优步开源的 Pyro 概率编程语言完成的，它是基于脸书的 PyTorch 深度学习库构建的。

**import** **random**

**import** **pandas** **as** **pd**
**import** **numpy** **as** **np**

**import** **torch**
**from** **pyro.infer.autoguide** **import** AutoDiagonalNormal 
**from** **pyro.infer** **import** SVI, Trace_ELBO, Predictive
**from** **pyro.distributions.distribution** **import** Distribution
**from** **torch.distributions** **import** constraints
**from** **pyro.distributions.torch_distribution** **import** TorchDistributionMixin, TorchDistribution
**import** **pyro.distributions** **as** **dist**
**import** **pyro**
**import** **seaborn** **as** **sns***# some utilities defined in https://github.com/agadetsky/pytorch-pl-variance-reduction/blob/master/utils.py*

**def** logcumsumexp(x, dim):
    *# slow implementation, but ok for now*
    **if** (dim != -1) **or** (dim != x.ndimension() - 1):
        x = x.transpose(dim, -1)

    out = []
    **for** i **in** range(1, x.size(-1) + 1):
        out.append(torch.logsumexp(x[..., :i], dim=-1, keepdim=**True**))
    out = torch.cat(out, dim=-1)

    **if** (dim != -1) **or** (dim != x.ndimension() - 1):
        out = out.transpose(-1, dim)
    **return** out

**def** reverse_logcumsumexp(x, dim):
    **return** torch.flip(logcumsumexp(torch.flip(x, dims=(dim, )), dim), dims=(dim, ))

**def** smart_perm(x, permutation):
    **assert** x.size() == permutation.size()
    **if** x.ndimension() == 1:
        ret = x[permutation]
    **elif** x.ndimension() == 2:
        d1, d2 = x.size()
        ret = x[
            torch.arange(d1).unsqueeze(1).repeat((1, d2)).flatten(),
            permutation.flatten()
        ].view(d1, d2)
    **elif** x.ndimension() == 3:
        d1, d2, d3 = x.size()
        ret = x[
            torch.arange(d1).unsqueeze(1).repeat((1, d2 * d3)).flatten(),
            torch.arange(d2).unsqueeze(1).repeat((1, d3)).flatten().unsqueeze(0).repeat((1, d1)).flatten(),
            permutation.flatten()
        ].view(d1, d2, d3)
    **else**:
        **ValueError**("Only 3 dimensions maximum")
    **return** ret**class** **PlackettLuce**(Distribution):
    *"""*
 *Plackett-Luce distribution*
 *"""*
    arg_constraints = {"logits": constraints.real}

    **def** __init__(self, logits):
        *# last dimension is for scores of plackett luce*
        super(PlackettLuce, self).__init__()
        self.logits = logits
        self.size = self.logits.size()

    **def** sample(self, sample_shape=torch.Size()):
        *# sample permutations using Gumbel-max trick to avoid cycles*
        **with** torch.no_grad():
            logits = self.logits.unsqueeze(0).expand(sample_shape, *self.size)
            u = torch.distributions.utils.clamp_probs(torch.rand_like(logits))
            z = self.logits - torch.log(-torch.log(u))
            samples = torch.sort(z, descending=**True**, dim=-1)[1]
        **return** samples

    **def** log_prob(self, samples):
        *# samples shape is: num_samples x self.size*
        *# samples is permutations not permutation matrices*
        **if** samples.ndimension() == self.logits.ndimension():  *# then we already expanded logits*
            logits = smart_perm(self.logits, samples)
        **elif** samples.ndimension() > self.logits.ndimension():  *# then we need to expand it here*
            logits = self.logits.unsqueeze(0).expand(*samples.size())
            logits = smart_perm(logits, samples)
        **else**:
            **raise** **ValueError**("Something wrong with dimensions")
        logp = (logits - reverse_logcumsumexp(logits, dim=-1)).sum(-1)
        **return** logp

好，让我们做一个小模拟。很简单。想象一下，我们有 5 名赛车手。我们将使用指数分布生成一些“价值”向量。我们看到，到目前为止，驾驶员 3 的评分最高，驾驶员 0 的评分最低。给定这些参数，然后我们使用 PL 分布生成“种族”的随机样本。我们这样做 10 次。

n_drivers = 5
torch.manual_seed(39)
theta = dist.Exponential(0.1).sample((n_drivers,))
ranks = list(zip(range(n_drivers),theta.numpy()))
[print(f'driver **{**i**}**: **{**j**:**0.3**}**') **for** i, j **in** sorted(ranks,key=**lambda** x: -x[1])]
pl = PlackettLuce(theta)
data = pl.sample(10)
print(data.numpy())# output
driver 3: 31.7
driver 4: 23.0
driver 2: 2.02
driver 1: 1.02
driver 0: 0.294
[[3 4 1 2 0]
 [3 4 1 2 0]
 [3 4 2 0 1]
 [3 4 2 1 0]
 [3 4 2 1 0]
 [3 4 1 0 2]
 [3 4 2 1 0]
 [3 4 1 2 0]
 [3 4 2 1 0]
 [3 4 0 2 1]]

你可以看到，在我们所有的 10 场模拟比赛中，3 号和 4 号车手分别获得了 1 号和 2 号。

现在，让我们做同样的模拟，但增加指数率参数，使驱动器更“均匀”。

n_drivers = 5
torch.manual_seed(39)
theta = dist.Exponential(1.0).sample((n_drivers,))
ranks = list(zip(range(n_drivers),theta.numpy()))
[print(f'driver **{**i**}**: **{**j**:**0.3**}**') **for** i, j **in** sorted(ranks,key=**lambda** x: -x[1])]
pl = PlackettLuce(theta)
data = pl.sample(10)
print(data.numpy())#output
driver 3: 3.17
driver 4: 2.3
driver 2: 0.202
driver 1: 0.102
driver 0: 0.0294
[[1 4 3 2 0]
 [3 1 4 2 0]
 [3 4 2 0 1]
 [4 3 2 0 1]
 [4 3 2 1 0]
 [3 4 0 1 2]
 [3 4 2 0 1]
 [3 4 1 2 0]
 [3 4 2 0 1]
 [4 0 3 2 1]]

你可以看到现在车手 4 赢得了一些比赛，甚至车手 1 也站在了领奖台上。太棒了。希望这开始对你有意义。我们基本上准备好了。我们只需要读入一些数据，创建模型，做一些推断。我们开始吧！

我从 17 个公开的模拟草稿中收集了数据。然而，并非所有的草案都模拟了两轮投票。这是可以的，因为即使在一个回合中仍然有有用的信息。我将很快解释我是如何处理这个问题的，但是让我们先看看数据。

df = pd.read_csv('data/mock-drafts.csv')
df.head()

作者图片

cols = df.columns[1:]
players_set = set()
**for** col **in** cols:
    **for** player **in** df[col].values:
        players_set.add(player)
players = {player:idx **for** idx,player **in** enumerate(players_set) **if** type(player) **is** str}
player_lut = {idx:player **for** player,idx **in** players.items() **if** type(player) **is** str}
print(player_lut)# output
{1: 'Kira Lewis', 2: 'Udoka Azubuike', 3: 'Tyrell Terry', 4: 'Saben Lee', 5: 'Aleksej Pokusevski', 6: 'Cassius Winston', 7: 'Austin Wiley', 8: 'Devon Dotson', 9: 'Naji Marshall', 10: 'Nick Richards', 11: 'Yam Madar', 12: 'Zeke Nnaji', 13: 'Kaleb Wesson', 14: 'Mason Jones', 15: 'Patrick Williams', 16: 'James Wiseman', 17: 'Tyrese Haliburton', 18: 'Theo Maledon', 19: 'LaMelo Ball', 20: 'Tyrique Jones', 21: 'Ty-Shon Alexander', 22: 'Vernon Carey', 23: 'Reggie Perry', 24: 'Xavier Tillman', 25: 'Filip Petrusev', 26: 'Leandro Bolmaro', 27: 'Josh Green', 28: 'Abdoulaye Ndoye', 29: 'Killian Hayes', 30: 'Tyler Bey', 31: 'Isaiah Stewart', 32: 'Lamine Diane', 33: 'Desmond Bane', 34: 'Lamar Stevens', 35: 'Aaron Nesmith', 36: 'Tyrese Maxey', 37: 'Anthony Edwards', 38: 'Yoeli Childs', 39: 'Omer Yurtseven', 40: 'Borisa Simanic', 41: 'Justinian Jessup', 42: 'Malachi Flynn', 43: 'Jahmius Ramsey', 44: 'Obi Toppin', 45: 'Markus Howard', 46: 'Tre Jones', 47: 'Robert Woodard', 48: 'Nate Hinton', 49: 'Isaac Okoro', 50: 'Cole Anthony', 51: 'Jaden McDaniels', 52: 'Jordan Nwora', 53: 'Deni Avdija', 54: 'Jalen Harris', 55: 'Rokas Jokubaitis', 56: 'Paul Reed', 57: 'Jay Scrubb', 58: 'Ashton Hagans', 59: 'Jalen Smith', 60: 'Kristian Doolittle', 61: 'Henri Drell', 62: 'CJ Elleby', 63: 'Isaiah Joe', 64: 'Paul Eboua', 65: 'Killian Tillie', 66: 'Kenyon Martin', 67: 'RJ Hampton', 68: 'Immanuel Quickley', 69: 'Sam Merrill', 70: 'Onyeka Okongwu', 71: 'Grant Riller', 72: 'Payton Pritchard', 73: 'Skylar Mays', 74: 'Karim Mane', 75: 'Giorgos Kalaitzakis', 76: 'Mamadi Diakite', 77: 'Nico Mannion', 78: 'Saddiq Bey', 79: 'Marko Simonovic', 80: 'Paul Ebouta', 81: 'Precious Achiuwa', 82: 'Cassius Stanley', 83: 'Daniel Oturu', 84: 'Elijah Hughes', 85: 'Josh Hall', 86: 'Devin Vassell'}

正如你所看到的，在这些模拟选秀中共有 86 名选手。现在我们将把模拟放入 Pyro 可以理解的数据格式中，一个包含 17 行的 Torch 张量。每一排将有所有 86 名球员在一些排列。现在我可以解释我是如何填充每个模拟的，这样无论一个模拟有 30 个还是 60 个玩家，每个都将被填充以包含所有的 86 个玩家。我假设一些非常简单的事情:如果一个模拟没有一个玩家，那么每个其他玩家可以被随机洗牌以填充整个模拟到位置 86。这是我能做的最简单的方法。更高级的方法是从这样的随机化开始，然后用基于模型推断评级的随机顺序来细化“反事实”填充。人们可以像这样不断迭代，直到达到某种平衡(基本上直到不再有改进)。

**import** **random**
n_players = len(player_lut)
list_data = []
**for** col **in** cols:
    mock = df[col]
    ranks = []
    players_cp = players.copy()
    **for** i **in** range(1,n_players+1):
        **try**:
            player = mock[i-1]
            players_cp.pop(player)
            ranks.append(players[player])
        **except** **KeyError**: 
            random_player = random.choice(list(players_cp.keys()))
            **if** players_cp.pop(random_player):
                ranks.append(players[random_player])
    list_data.append(ranks)
data = torch.tensor(list_data)-1
print(data)# output
tensor([[18, 15, 36,  ..., 44, 73,  6],
        [36, 52, 15,  ...,  8, 53, 78],
        [18, 16, 36,  ..., 75, 61, 63],
        ...,
        [18, 15, 36,  ..., 31,  3,  6],
        [36, 18, 69,  ..., 72, 51, 63],
        [36, 15, 18,  ..., 37, 44, 19]])

让我们看看其中的一行:

[player_lut[i+1] **for** i **in** data.numpy()[0]][0:10]# output
['LaMelo Ball',
 'James Wiseman',
 'Anthony Edwards',
 'Deni Avdija',
 'Obi Toppin',
 'Tyrese Haliburton',
 'Isaac Okoro',
 'Devin Vassell',
 'Onyeka Okongwu',
 'Killian Hayes']

好了，现在是激动人心的部分。让我们定义模型:

n_players = 86

**def** model(data):
    sigma = pyro.sample("sigma", dist.Uniform(0,1))
    theta = pyro.sample("theta", dist.Cauchy(torch.zeros(n_players,),sigma*torch.ones(n_players,)))
    **with** pyro.plate("data", data.shape[-2]):
        pyro.sample("obs", PlackettLuce(theta), obs=data)

真的就是这样。我使用柯西分布，因为如果你考虑 NBA 人才的分布，它似乎是“厚尾”的。换句话说，我们希望我们的模型能够处理有“超级明星”或一致同意的最佳人选的情况。如果我们让我们的分布过于平坦，这将更像一场赛马，任何一匹马都有机会第一个到达终点(或多或少)。现在我们进行推理步骤，根据 Pyro 的说法建立一个向导和优化器来进行随机变分推理(在未来，我计划研究更精确的 MCMC 采样方法，但是现在应该可以了)。

guide = AutoDiagonalNormal(model)
adam = pyro.optim.Adam({"lr": 0.03})
svi = SVI(model, guide, adam, loss=Trace_ELBO())
pyro.clear_param_store()
**for** j **in** range(5000):
    *# calculate the loss and take a gradient step*
    loss = svi.step(data)
    **if** j % 100 == 0:
        print("[iteration **%04d**] loss: **%.4f**" % (j + 1, loss / data.shape[-2]))

这是我们元模拟的结果！

1\. Anthony Edwards (12.2)
2\. LaMelo Ball (12.0)
3\. Deni Avdija (10.3)
4\. James Wiseman (10.3)
5\. Obi Toppin (9.9)
6\. Onyeka Okongwu (8.0)
7\. Killian Hayes (7.9)
8\. Tyrese Haliburton (7.6)
9\. Isaac Okoro (7.6)
10\. Devin Vassell (6.7)
11\. Patrick Williams (6.0)
12\. Aaron Nesmith (5.2)
13\. Saddiq Bey (4.7)
14\. Kira Lewis (4.6)
15\. Precious Achiuwa (4.5)
16\. RJ Hampton (4.1)
17\. Jalen Smith (4.0)
18\. Cole Anthony (3.9)
19\. Tyrese Maxey (3.7)
20\. Josh Green (3.0)
21\. Aleksej Pokusevski (2.6)
22\. Leandro Bolmaro (2.4)
23\. Isaiah Stewart (2.3)
24\. Jaden McDaniels (2.2)
25\. Theo Maledon (2.0)
26\. Jahmius Ramsey (1.7)
27\. Tyrell Terry (1.6)
28\. Desmond Bane (1.3)
29\. Nico Mannion (1.2)
30\. Zeke Nnaji (1.1)
31\. Malachi Flynn (1.0)
32\. Tre Jones (0.9)
33\. Devon Dotson (0.8)
34\. Paul Reed (0.5)
35\. Vernon Carey (0.5)
36\. Isaiah Joe (0.5)
37\. Elijah Hughes (0.5)
38\. Cassius Stanley (0.5)
39\. Daniel Oturu (0.4)
40\. Killian Tillie (0.4)
41\. Cassius Winston (0.4)
42\. Grant Riller (0.3)
43\. Udoka Azubuike (0.3)
44\. Tyler Bey (0.3)
45\. Reggie Perry (0.2)
46\. Xavier Tillman (0.2)
47\. Payton Pritchard (0.1)
48\. Jay Scrubb (0.1)
49\. Abdoulaye Ndoye (0.1)
50\. Immanuel Quickley (-0.0)
51\. Ty-Shon Alexander (-0.0)
52\. Ashton Hagans (-0.0)
53\. Jordan Nwora (-0.1)
54\. Robert Woodard (-0.1)
55\. Nick Richards (-0.1)
56\. Yam Madar (-0.1)
57\. Skylar Mays (-0.2)
58\. Nate Hinton (-0.2)
59\. Paul Ebouta (-0.2)
60\. Filip Petrusev (-0.3)
61\. Justinian Jessup (-0.3)
62\. Lamar Stevens (-0.3)
63\. Marko Simonovic (-0.3)
64\. Sam Merrill (-0.3)
65\. Josh Hall (-0.4)
66\. Karim Mane (-0.4)
67\. Kaleb Wesson (-0.4)
68\. Kristian Doolittle (-0.5)
69\. Tyrique Jones (-0.5)
70\. Borisa Simanic (-0.5)
71\. Markus Howard (-0.5)
72\. Kenyon Martin (-0.5)
73\. CJ Elleby (-0.6)
74\. Rokas Jokubaitis (-0.6)
75\. Lamine Diane (-0.6)
76\. Yoeli Childs (-0.6)
77\. Jalen Harris (-0.6)
78\. Mason Jones (-0.7)
79\. Paul Eboua (-0.7)
80\. Saben Lee (-0.7)
81\. Mamadi Diakite (-0.7)
82\. Omer Yurtseven (-0.8)
83\. Henri Drell (-0.8)
84\. Austin Wiley (-0.8)
85\. Naji Marshall (-0.8)
86\. Giorgos Kalaitzakis (-0.9)

好吧，还不错，但这还不是最酷的部分。最酷的部分是现在我们有了每个玩家的概率分布，所以我们可以去模拟成千上万的模拟选秀。我们开始吧！基本上，我们只是对后验分布进行采样(来自 Pyro guide 函数),并生成一个样本模拟汇票订单。然后我们做那 10，000(！)次，并将它们全部放入熊猫数据框中。

rank_hist = []
**for** m **in** range(10_000):
    theta = guide.forward()['theta']
    pl = PlackettLuce(theta)
    sample = pl.sample(1).numpy()[0]
    **for** i, j **in** enumerate(sample.squeeze()):
        **if** i<=59:
            rank_hist.append((i+1,player_lut[j+1]))
        **else**:
            rank_hist.append((61,player_lut[j+1]))
rank_df = pd.DataFrame.from_records(rank_hist,columns=['rank','player'])

所以现在我们有了每个球员的模拟选秀位置的分布(直方图)。例如，让我们看看安东尼·爱德华兹是在哪里被预言的:

sns.set_style("ticks", {"xtick.minor.size": 1, "ytick.major.size": 8})
players = ['Anthony Edwards']
sns.displot(data=rank_df[rank_df.player.isin(players)],x='rank',kind='hist',hue='player')

作者图片

理解这种分布的一个更简单的方法是看累积概率。还记得我说过这个模型可以以类似生存分析的方式使用吗？这就是我所说的:

sns.set_style("ticks", {"xtick.minor.size": 1, "ytick.major.size": 8})
players = ['Anthony Edwards']
sns.displot(data=rank_df[rank_df.player.isin(players)],x='rank',kind='ecdf',hue='player')

作者图片

从这里我们看到爱德华兹有大约 45%的机会被选为第一名，大约 80%的机会进入前两名，以此类推。相当整洁！现在让我们来看看一个预测范围更广的人，Saddiq Bey:

sns.set_style("ticks", {"xtick.minor.size": 1, "ytick.major.size": 8})
players = ['Saddiq Bey']
sns.displot(data=rank_df[rank_df.player.isin(players)],x='rank',kind='ecdf',hue='player')

作者图片

哦，那不是很酷吗？确实范围很广。让我们来看看一个更加边缘的球员，尼克·理查兹:

sns.set_style("ticks", {"xtick.minor.size": 1, "ytick.major.size": 8})
players = ['Nick Richards']
sns.displot(data=rank_df[rank_df.player.isin(players)],x='rank',kind='ecdf',hue='player')

作者图片

好的，所以我应该解释一下，在模拟中，如果一个球员没有被选进前 60，我就把他放在 61，这样他基本上就是一个未被选中的自由球员(UDFA)。这些图有点乱，但在我完成之前，这里有一个图表中有前 10 名的人。

sns.set_style("ticks", {"xtick.minor.size": 1, "ytick.major.size": 8})
players = ['Anthony Edwards', 'LaMelo Ball', 'Deni Avdija', 'James Wiseman', 'Obi Toppin', 'Onyeka Okongwu',
          'Killian Hayes', 'Tyrese Haliburton', 'Isaac Okoro', 'Devin Vassell']
sns.displot(data=rank_df[rank_df.player.isin(players)],x='rank',kind='ecdf',hue='player')

作者图片

嗯，我现在就知道这么多了。随着选秀临近(11 月 18 日)，请(通过我的推特)关注更多与此相关的事情。如果你有任何建议(或修正)让我知道！以下是该项目的 github repo 的链接:

预测分析:Python 中的贝叶斯线性回归

原文：https://towardsdatascience.com/bayesian-linear-regression-project-in-python-forecast-water-consumption-under-the-impact-of-cea62c2693e4?source=collection_archive---------18-----------------------

Alex Samuels 在 Unsplash 上的照片

数据科学项目的贝叶斯线性回归基础概述

在这篇文章中，我想更多地关注贝叶斯线性回归理论，并用 Python 实现一个数据科学项目的建模。整个项目是关于预测未来三十年气候变化影响下的城市用水量。完整版本的代码可以在 GitHub 的 Jupyter 笔记本上找到，我鼓励你去看看。此外，研究结果发表在可持续发展的城市和社会上。

为什么是贝叶斯统计？

统计学有两个主要分支，描述统计和推断统计。描述性统计提供数据的简明摘要，而推断性统计发现模式并对数据做出推断。但是为什么我要谈论这两个术语呢？😀嗯，在推断统计学中有两个对我们很重要的一般“哲学”;频率主义和贝叶斯主义。

频率统计和贝叶斯统计通过两种不同的方法来划分统计世界，以定义概率。在频率统计中，概率与重复事件的频率有关，而在贝叶斯统计中，概率与我们自身的确定性和不确定性有关。因此，结果是，频率主义者认为模型是固定的，数据围绕它们变化，但贝叶斯主义者谈论观察到的数据是固定的，模型可以围绕它们变化。

所以这个问题的答案是“为什么是贝叶斯统计？”是，有了贝叶斯统计，我们可以用概率来表示任何事件或假设中的不确定性。

贝叶斯定理

贝叶斯定理是贝叶斯统计中的关键规则[1]。如果 H 是假设，D 是观测数据，贝叶斯定理如下

贝叶斯定理

其中 P(H|D) 表示假设 D 发生的条件下，假设的后验概率。 P(H) 是假设的先验概率， P(D) 是观测数据的边际(总)概率，实际上是一个归一化常数， P(D|H) 是给定假设的数据的概率(可能性)。

资料组

加拿大魁北克省的 Brossard 市被选为研究地点。我获得了 2011 年 1 月至 2015 年 10 月期间的日城市用水量、 q 的时间序列。对于同一时间段，我还得到了气候变量的时间序列:日最低气温 θ ，日最高气温 t ，总降水量 p 。这些气候变量的测量是在皮埃尔·艾略特·特鲁多国际机场(YUL)的加拿大环境站进行的。

贝叶斯线性回归

贝叶斯统计是一种强大的概率建模技术，已被广泛用于统计建模，包括线性回归模型，以预测系统[2，3，4，5]。线性回归模型表示为

线性回归模型

其中 q 为因变量(目标变量) β_0 为截距； β_1 、 β_2 、 β_3 为权重(称为模型参数)； θ 、 t 和 p 称为预测变量；并且 ε 是表示随机噪声或模型方程中不包括的变量的影响的误差项。该等式可以改写为

其中β=(β_ 0， β_1 ， β_2 ，β_ 3)；以及X*=*(1，q ， Q ，p)【2，6】。线性回归是使用概率分布而不是点估计来预测 q 的。因此， q 不是作为单个值来估计的，而是假设从概率分布中提取的。响应从正态(高斯)分布 N 中取样的贝叶斯线性回归模型为

输出 q 由以平均值和方差为特征的正态分布产生。正态分布的平均值是回归系数矩阵()乘以预测值矩阵( X )。方差是标准差的平方， σ 。

贝叶斯线性回归模型提供了预测变量中不确定性的表示，并确定了模型参数的后验分布。不仅响应是由概率分布生成的，而且模型参数也被假定来自分布。模型参数的后验概率来自贝叶斯定理:

模型参数的后验概率

其中P*(β*| q，* X )是给定预测值和自变量的模型参数的后验概率分布； P(q| β，X ) 是数据的可能性；*P(β|X)是参数的先验概率，方程的分母是根据全概率定律可以找到的边际概率。

马尔可夫链蒙特卡罗

在实践中，很难计算连续值的边际可能性；计算精确的后验分布是困难的。作为一种解决方案，一种采样方法马尔可夫链蒙特卡罗被用于近似后验概率，而无需计算边际似然[3，5，7]。蒙特卡洛是一种抽取随机样本的通用技术，马尔可夫链意味着抽取的下一个样本仅基于上一个样本值。通过引入更多样本，后验概率的近似值将最终收敛于 β_1 、 β_2 和 β_3 的实际后验概率分布。

作为应用马尔可夫链蒙特卡罗的起点，定义覆盖 β_1 、 β_2 和 β_3 所有可能值的参数空间。然后，将每个参数的先验概率定义为正态分布。接下来，计算每个可能参数的可能性。最后，计算参数空间中任何给定点的先验×似然。

使用不掉头算法【8】实现马尔可夫链蒙特卡罗。当所涉及的变量是连续的并且不需要设置步数时，这种算法是有效的。这是优于哈密尔顿蒙特卡罗的优点，哈密尔顿蒙特卡罗采用由一阶梯度信息通知的一系列步骤，并且对步骤的数量敏感。

生成用于被接受的移动的一组参数值(如果建议的移动未被接受，则重复先前的值),并且在多次重复之后，找到分布的经验近似。最终，执行贝叶斯线性回归的结果是基于数据和先验的可能模型参数的概率密度函数。

模型性能测量

整个数据集被分成训练集和测试集。训练集包含 75%的数据，用于构建模型；而测试集包含 25%，用于验证预测的准确性。此外，两个拟合优度指标用于评估本研究中开发的模型的性能。这些度量是平均绝对误差和均方根误差。平均绝对误差是预测值和实际值之差的绝对值的平均值，均方根误差是预测值和实际值之差的平均值的平方根。

均方根误差

绝对平均误差

Python 代码

正如我之前提到的，我在这篇文章中的重点是贝叶斯线性回归模型。因此，我跳过数据预处理这一步。贝叶斯线性回归模型的输入数据集如下

输入数据集

• UWC:城市用水量
• Max_T:最高温度
• 最低温度
• 总降水量

⚠️在创建模型之前，请确保您已经处理了缺失值，因为贝叶斯线性回归模型不接受带有缺失点的数据。

在 PyMC3 中创建贝叶斯线性回归模型

首先，我使用 sklearn 库将预处理数据集(df)拆分为 75%训练和 25%测试。

注: X_train 和 X_test 包含目标变量 UWC。

*from sklearn.model_selection import train_test_splitX_train, X_test, y_train, y_test = train_test_split(df, df['UWC'],
test_size = 0.25,random_state = 42)*

在这个项目中，我在 PyMC3 中创建了贝叶斯线性回归模型。让马尔可夫链蒙特卡罗算法从后验概率中抽取样本来近似每个模型参数的后验概率。

为了创建模型，我将线性回归公式、参数的先验族和数据传递给 GLM 。为了执行马尔可夫链蒙特卡罗抽样，我指定样本数为 2000，链数为 2，调优步数为 500。

*import pymc3 as pm
from pymc3 import traceplot**# Formula for Bayesian Linear Regression** 
formula = ‘UWC ~ ‘ + ‘ + ‘.join([‘%s’ % variable for variable in X_train.columns[1:]])
print(formula)**# Context for the model**
with pm.Model() as normal_model:

    **# The prior for the model parameters will be a normal distribution**
    family = pm.glm.families.Normal()**# Creating the model requires a formula and data (and optionally a family)**
    pm.GLM.from_formula(formula, data = X_train, family = family)**# Perform Markov Chain Monte Carlo sampling**
    normal_trace = pm.sample(draws=2000, chains = 2, tune = 500)*

所以，UWC 是最高温度、最低温度和总降水量的函数。

将后验分布绘制成直方图

*pm.plot_posterior(normal_trace)*

(a)截距β_0 的后验直方图；(b)对于 θ 的回归系数β_ 1*；©回归系数 β_2 为t；(d)对于 p 的回归系数β_ 3；标准偏差 σ 。*

*结果显示为 β_0 、 β_1 、 β_2 、 β 、 _ 、 *3、和 σ 的后验直方图。直方图显示 95%的最高后验密度(HPD)，这是参数的可信区间。贝叶斯统计中的可信区间相当于置信区间。比如 0.034<σ0.037 的概率是 95%。

计算梅和 RMSE

为了计算 MAE 和 RMSE，我取 UWC 概率分布的中值，并将其作为对 evaluate_prediction 函数的预测。

***# Define a function to calculate MAE and RMSE**
def evaluate_prediction(prediction, true):
    mae = np.mean(abs(predictions - true))
    rmse = np.sqrt(np.mean((predictions - true) ** 2))

    return mae, rmsemedian_pred = X_train['UWC'].median()
median_preds = [median_pred for _ in range(len(X_test))]
true = X_test['UWC']**# Display mae and rmse**
mae, rmse = evaluate_prediction(median_preds, true)
print('Mean Absolute Error: {:.4f}'.format(mae))
print('Root Mean Square Error: {:.4f}'.format(rmse))*

对模型性能的评估表明，预测模型具有良好的精度。

结论

这篇文章是一个数据科学项目的贝叶斯线性回归的简要概述。我希望它能帮助你更好地理解贝叶斯线性回归的基本原理。😊

非常感谢您的反馈。你可以在 LinkedIn 上找到我。

参考

*[1]赫克曼博士、盖格博士和奇克林博士(1995 年)。学习贝叶斯网络:知识和统计数据的结合。机器学习， 20 (3)，197–243。【https://doi.org/10.1023/A:1022623210503 *

[2]韩，贾伟和裴，简和坎伯，M. (2011)。数据挖掘:概念和技术。检索自http://myweb . saban ciuniv . edu/rdehkharghani/files/2016/02/The-mor gan-Kaufmann-Series-in-Data-Management-Systems-Jia Wei-Han-Micheline-Kamber-Jian-Pei-Data-Mining。-概念和技术-第三版-摩根-考夫曼-2011.pdf

[3] Mudgal，a .，Hallmark，s .，Carriquiry，a .，和 Gkritza，K. (2014 年)。环形交叉路口的驾驶行为:分层贝叶斯回归分析。交通研究 D 部分:交通与环境， 26 ，20–26。https://doi.org/10.1016/j.trd.2013.10.003

[4]a . e .拉夫特里、d .马迪根和 J. A .赫廷(1997 年)。线性回归模型的贝叶斯模型平均。美国统计协会杂志， 92 (437)，179–191。检索自https://www . tandfonline . com/doi/ABS/10.1080/01621459 . 1997 . 10473615

[5]袁晓春，孙晓霞，赵，魏，米，张，王，魏延明(2017)。2030 年中国区域能源需求预测:贝叶斯方法。资源、保护和回收，127(5 月)，85–95。https://doi.org/10.1016/j.resconrec.2017.08.016

[6]t .塔勒布和 m .卡杜尔(2017 年)。无线传感器网络中提高能量效率的分层凝聚分簇方案。交通与电信， 18 (2)，128–138。https://doi.org/10.1515/ttj-2017-0012

[7] Godsill，S. J. (2001 年)。模型不确定性的马尔可夫链蒙特卡罗方法之间的关系。计算与图形统计杂志， 10 (2)，230–248。https://doi.org/10.1198/10618600152627924

[8]医学博士霍夫曼和盖尔曼(2011 年)。不掉头采样器:在哈密顿蒙特卡罗中自适应地设置路径长度。(2008), 1–30.从 http://arxiv.org/abs/1111.4246取回

使用 PyMC3 的贝叶斯逻辑回归

原文：https://towardsdatascience.com/bayesian-logistic-regression-with-pymc3-8e17c576f31a?source=collection_archive---------11-----------------------

这是一系列文章中的另一篇(参见此处和此处的其他相关文章),一般来说是关于概率编程的，特别是关于 PyMC3 的。在我们之前的文章中，我们解释了 PyMC3 如何帮助统计推断。在本文中，我们将使用 PyMC3 解决一个端到端的分类问题。更准确地说，我们将使用 PyMC3 通过以下公共数据集进行贝叶斯逻辑回归:

https://archive . ics . UCI . edu/ml/datasets/Occupancy+Detection+

该数据集包含几个变量，如光线、温度、湿度和二氧化碳水平。目标是从这些变量中检测房间的占用情况。

首先，我们需要加载几个相关的包。

**import** **arviz** **as** **az**
**import** **matplotlib.pyplot** **as** **plt**
**import** **numpy** **as** **np**
**import** **pandas** **as** **pd**
**import** **pymc3** **as** **pm**
**import** **seaborn**
**import** **theano.tensor** **as** **tt**
**import** **warnings**
**from** **IPython.core.pylabtools** **import** figsize
**import** **seaborn** **as** **sns**
**from** **sklearn.metrics** **import** (roc_curve, roc_auc_score, confusion_matrix, accuracy_score, f1_score, 
                             precision_recall_curve) 
**from** **sklearn.metrics** **import** confusion_matrix

接下来，我们加载数据集。

df=pd.read_csv('datatest.txt')
df.sample(5)

为了方便起见，我们将日期变量转换为日期时间对象。

df['date']=pd.to_datetime(df['date'])

第 1 部分:探索性数据分析

首先，让我们看一下数据集的概况。

df.describe()

数据集概述

我们看到这个数据集中有 2655 个样本。此外，没有丢失值。让我们也看看这个数据集的时间框架。

df['date'].describe()count                    2665
unique                   2665
top       2015-02-03 07:25:59
freq                        1
first     2015-02-02 14:19:00
last      2015-02-04 10:43:00
Name: date, dtype: object

所以我们的数据只在三天内就收集好了。接下来，我们将探讨我们的变量及其关系。首先，让我们画出温度变量。

figsize(12.5, 4)
plt.hist(df['Temperature'], bins=40, density=**True**, label='Temperature')
plt.xlabel('Temperature')
plt.title('Distribution of temperature')
plt.show()

该图显示温度具有重尾分布。湿度变量呢？

figsize(12.5, 4)
plt.hist(df['Humidity'], bins=50, density=**True**, label='Humidity')
plt.xlabel('Humidity')
plt.title('Distribution of Humidity')
plt.show()

有趣的是，在 22.5 和 25 附近有两个峰值。我们还对不同日子里的光照变化感兴趣。

figsize(12.5, 4)
sns.boxplot(x=df['date'].dt.day,y=df['Light'], orient='v')
plt.xlabel('Day')
plt.title('Boxplot for Light during different days')
plt.show()

我们看到光的分布在这三天中几乎是相同的。接下来，让我们看看二氧化碳水平。

figsize(12.5, 4)
sns.boxplot(x=df['date'].dt.day,y=df['CO2'], orient='v')
plt.xlabel('Day')
plt.title('Boxplot for CO2 level during different days')
plt.show()

这些分布明显不同。2015 年 2 月 4 日有很多异常值。最后，我们将深入研究湿度比变量。

figsize(12.5, 4)
sns.boxplot(x=df['date'].dt.day,y=df['HumidityRatio'], orient='v')
plt.xlabel('Day')
plt.title('Boxplot for Humidity Ratio level during different days')
plt.show()

这看起来非常类似于二氧化碳水平的箱线图。也许，二氧化碳水平和湿度之间有很强的相关性。我们可以用这两个变量的散点图来验证。

ax=sns.scatterplot(df['CO2'], df['HumidityRatio'], style=df['date'].dt.day)

二氧化碳水平和湿度比的散点图

事实上，二氧化碳水平和湿度比之间有很强的线性关系。让我们看一下变量之间的关系。这可以通过 seaborn 的 pair plot 函数来完成。

ax=seaborn.pairplot(df)

数据集中所有要素的配对图

这直观地显示了在以下对之间存在强线性关系:CO2 和温度、CO2 和湿度、湿度和湿度比、湿度比和 CO2。我们甚至可以通过绘制热图来量化这些关系。

corr=df.iloc[:, 1:-1].corr()
*# Generate a mask for the upper triangle*
mask = np.triu(np.ones_like(corr, dtype=np.bool))
*# Set up the matplotlib figure*
f, ax = plt.subplots(figsize=(11, 9))*# Draw the heatmap with the mask and correct aspect ratio*
ax=sns.heatmap(corr, mask=mask,
            square=**True**, linewidths=.5, cbar_kws={"shrink": .5})

清楚地表明，两对湿度-湿度比和湿度比-CO2 表现出最强的线性关系。

基于 PyMC3 的贝叶斯逻辑模型

我们将建立几个机器学习模型，根据其他变量对入住情况进行分类。

回想一下，我们有一个二元决策问题。换句话说，我们的目标变量假设遵循伯努利随机变量，p 由下式给出:

其中，var 是我们在模型中使用的所有变量的集合，logit 是逻辑函数。

为了建立贝叶斯逻辑回归模型，我们首先要对每个参数进行先验分布。这些先验的选择将影响结果(尽管有更多的数据，它们可能会“收敛”到相同的分布。)

一旦我们的先验被指定，PyMC3 将使用马尔可夫链蒙特卡罗模拟及其推广在数值上近似后验分布。然后，我们可以使用这些后验样本进行推断。

先验分布均匀的第一个模型

由于我们没有关于这些参数的先验知识，我们可以假设它们可以是任何东西。换句话说，我们假设所有的β_var 服从一个上下界都很大的均匀分布。为了捕捉一张大网，我们对均匀分布使用大的上下界。

 lower=-10**6
higher=10**6
**with** pm.Model() **as** first_model:
    *#priors on parameters*
    beta_0=pm.Uniform('beta_0', lower=lower, upper= higher)
    beta_temp=pm.Uniform('beta_temp', lower, higher)
    beta_humid=pm.Uniform('beta_humid', lower, higher)
    beta_light=pm.Uniform('beta_light', lower, higher)
    beta_co2=pm.Uniform('beta_co2', lower, higher)
    beta_humid_ratio=pm.Uniform('beta_humid_ration', lower, higher)

    *#the probability of belonging to class 1*
    p = pm.Deterministic('p', pm.math.sigmoid(beta_0+beta_temp*df['Temperature']+
                               beta_humid*df['Humidity']+
                               beta_light*df['Light']+
                               beta_co2*df['CO2']+
                               beta_humid_ratio*df['HumidityRatio']))
**with** first_model:
    *#fit the data* 
    observed=pm.Bernoulli("occupancy", p, observed=df['Occupancy'])
    start=pm.find_MAP()
    step=pm.Metropolis()

    *#samples from posterior distribution* 
    trace=pm.sample(25000, step=step, start=start)
    burned_trace=trace[15000:]

这可能需要一段时间来运行。一旦完成，我们就可以绘制样本。

pm.traceplot(burned_trace)
plt.show()

我们得出结论，我们的算法确实收敛。我们可以计算这些后验分布的平均值。

coeffs=['beta_0', 'beta_temp', 'beta_humid', 'beta_light', 'beta_co2', 'beta_humid_ration']
d=dict()
**for** item **in** coeffs:
    d[item]=[burned_trace[item].mean()]

result_coeffs=pd.DataFrame.from_dict(d)    
result_coeffs
*#coeff_result=pd.DataFrame(d)* 
*#coeff_result*

与频率统计相比，贝叶斯统计的一个优势是我们拥有的不仅仅是平均值。特别是，我们可以计算这些参数的 95% 高密度区间%2C%202015-,4.3.,highest%20density%20interval%2C%20abbreviated%20HDI.&text=Moreover%2C%20the%20probability%20density%20of,any%20x%20outside%20those%20limits.)。

pm.stats.hpd(burned_trace['beta_0'])
coeffs=['beta_0', 'beta_temp', 'beta_humid', 'beta_light', 'beta_co2', 'beta_humid_ration']
interval=dict()**for** item **in** coeffs:
    interval[item]=pm.stats.hpd(burned_trace[item]) *#compute 95% high density interval*

result_coeffs=pd.DataFrame.from_dict(interval).rename(index={0: 'lower', 1: 'upper'})
result_coeffs

注意，湿度比的系数明显大于其他系数。这并不一定意味着这个变量更重要。如果我们仔细观察数据，我们会发现这个变量取的值非常小。

此外，我们可以探索不同参数之间的关系。例如，让我们来看看β_ CO2 和β_ wet _ ratio 系数。

figsize(12.5, 12.5)
seaborn.jointplot(burned_trace['beta_co2'], burned_trace['beta_humid_ration'], kind="hex") *#color="#4CB391")*
plt.xlabel("beta_co2")
plt.ylabel("beta_humid_ratio");

CO2 和湿度比系数的散点图

曲线图显示这两个系数负相关。注意，CO2 水平和湿度比是正相关的。

回想一下，在经典的逻辑回归中，我们通过最大后验估计 (MAP 解)来寻找最佳参数。换句话说，最佳拟合参数由下式给出

其中，𝑝(𝜃|𝐷)p(θ|D 是给定数据时θ的后验分布，p(D|θ)是似然函数，p(θ)是θ的先验分布。

请注意，由于我们在第一个模型中使用均匀分布，我们可以预期我们的 MAP 解决方案应该与 MLE 解决方案(最大似然估计)一致，后者对应于 frequentist logistic 回归。我们可以使用 Scikit-Learn 库来测试这个语句。首先，我们使用 MAP 计算系数。

coeffs=['beta_0', 'beta_temp', 'beta_humid', 'beta_light', 'beta_co2', 'beta_humid_ration']
d=dict()
**for** item **in** coeffs:
    d[item]=[float(start[item])]

map_coeffs=pd.DataFrame.from_dict(d)    
map_coeffs

接下来，我们使用经典的逻辑回归计算β系数。

**from** **sklearn.linear_model** **import** LogisticRegressionX=df.iloc[:, 1: -1]
y=df['Occupancy']
logit=LogisticRegression()
logit_model=logit.fit(X,y)
pd.DataFrame(logit_model.coef_, columns=X.columns)

是啊！两种方法的系数几乎相同。

接下来，我们讨论了我们的模型的预测能力，并与经典的逻辑回归进行了比较。我们用经典方法记录预测。

logit_prediction=logit_model.predict(X)

为了使用贝叶斯逻辑模型进行预测，我们通过对样本值进行平均来计算𝑦_score。

*#compute the average probability of predicting 1\.* 
y_score = np.mean(burned_trace['p'], axis=0)*#histogram of the distribution*
figsize(12.5,4)
plt.hist(y_score, bins=40, density=**True**)
plt.xlabel('Probability')
plt.ylabel('Frequency')
plt.title('Distribution of $y_score$')
plt.show()

y 分数的分布

有趣的是，p 的大部分集中在 0 和 1 附近。我们也可以使用 y_score 进行预测。

first_model_prediction=[1 **if** x >0.5 **else** 0 **for** x **in** y_score]

让我们通过计算混淆矩阵来评估我们的模型的性能。

first_model_confussion_matrix =confusion_matrix(df['Occupancy'], first_model_prediction)
first_model_confussion_matrixarray([[1639,   54],
       [   3,  969]])

这很好。我们甚至还可以通过其他指标来量化性能。

**import** **sklearn**
**from** **sklearn.metrics** **import** classification_reportprint(sklearn.metrics.classification_report(y, first_model_prediction))precision    recall  f1-score   support 0       1.00      0.97      0.98      1693
           1       0.95      1.00      0.97       972 accuracy                           0.98      2665
   macro avg       0.97      0.98      0.98      2665
weighted avg       0.98      0.98      0.98      2665

我们也可以计算曲线下的面积。

pred_scores = dict(y_true=df['Occupancy'],y_score=y_score)
roc_auc_score(**pred_scores)0.99358530283253

所以，我们的模型表现得相当好。我们来和经典的 logistic 回归比较一下。

print(sklearn.metrics.classification_report(y, logit_prediction))precision    recall  f1-score   support 0       1.00      0.97      0.98      1693
           1       0.95      1.00      0.97       972 accuracy                           0.98      2665
   macro avg       0.97      0.98      0.98      2665
weighted avg       0.98      0.98      0.98      2665

他们是一样的！然而，使用贝叶斯模型，我们获得了更多的信息，因此我们对自己的估计更有信心。

先验分布为正态分布的不同模型。

现在，让我们使用一组不同的先验知识来训练我们的模型。例如，我们可以假设系数遵循正态分布。

**with** pm.Model() **as** second_model: 
    #priors with normal distribution beta_0=pm.Normal('beta_0', mu=0, sd=10**4)
    beta_temp=pm.Normal('beta_temp', mu=0, sd=10**4)
    beta_humid=pm.Normal('beta_humid', mu=0, sd=10**4)
    beta_light=pm.Normal('beta_light', mu=0, sd=10**4)
    beta_co2=pm.Normal('beta_co2', mu=0, sd=10**4)
    beta_humid_ratio=pm.Normal('beta_humid_ration', mu=0, sd=10**4) #probability of belonging to class 1 p = pm.Deterministic('p', pm.math.sigmoid(beta_0+beta_temp*df['Temperature']+
                               beta_humid*df['Humidity']+
                               beta_light*df['Light']+
                               beta_co2*df['CO2']+
                               beta_humid_ratio*df['HumidityRatio']))#fit observed data into the model**with** second_model:
    observed=pm.Bernoulli("occupancy", p, observed=df['Occupancy'])
    start=pm.find_MAP()
    step=pm.Metropolis()
    second_trace=pm.sample(25000, step=step, start=start)
    second_burned_trace=second_trace[15000:]
pm.traceplot(second_burned_trace)
plt.show()

我们再次看到算法确实收敛了。让我们计算 MAP 解的β系数。

coeffs=['beta_0', 'beta_temp', 'beta_humid', 'beta_light', 'beta_co2', 'beta_humid_ration']
d=dict()
**for** item **in** coeffs:
    d[item]=[float(start[item])]

second_map_coeffs=pd.DataFrame.from_dict(d)    
second_map_coeffs

它们与我们在第一个模型中得到的非常接近。更进一步，因为我们在贝叶斯框架中，我们甚至可以比较两个模型的后验分布。例如，让我们看看截距变量。

figsize(12.5,4)
plt.hist(burned_trace['beta_0']-second_burned_trace['beta_0'], bins=40, density=**True**)
plt.title('Distribution of the difference between beta_0')
plt.legend()
plt.show()

虽然 MAP 解对β_0 给出了相同的估计，但我们看到两个后验概率相当不同。让我们也比较一下这两个模型之间β_temp 的后验分布。

figsize(12.5,4)
plt.hist(burned_trace['beta_temp'], label='First model', bins=40, density=**True**)
plt.hist(second_burned_trace['beta_temp'], bins=40, label='Second model', density=**True**)
plt.title('Distribution of of beta_temp')
plt.legend()
plt.show()

实际上差别很小。接下来，让我们计算第二个模型的预测能力。

second_y_score = np.mean(second_burned_trace['p'], axis=0)
second_model_prediction=[1 **if** x >0.5 **else** 0 **for** x **in** second_y_score]
second_model_confussion_matrix =confusion_matrix(df['Occupancy'], second_model_prediction)
second_model_confussion_matrixarray([[1639,   54],
       [   3,  969]])

这与我们从第一个模型中得到的结果相同。我们可以检查一下 y_score 和 second_y_score 几乎是一样的。

figsize(12.5,4)
plt.hist(y_score-second_y_score, bins=40)
plt.title('Distribution of the difference between y_score and second_y_score')
plt.ylabel('Frequency')
plt.show()

内置广义线性模型(GLM)的模型。

在前面的章节中，我们使用实际操作的方法来构建我们的模型。这很容易，因为我们只有几个变量。当变量数量很大的时候，就不会很实用了。幸运的是，PyMC3 有一个内置的广义线性模型，在这个模型中，一切都将自动化。让我们使用这个内置模型来拟合我们的数据。

**with** pm.Model() **as** third_model:
    pm.glm.GLM.from_formula('Occupancy ~ Temperature + Humidity + Light + CO2 + HumidityRatio',
                            df,
                            family=pm.glm.families.Binomial())
    third_trace = pm.sample(25000, tune=10000, init='adapt_diag')pm.traceplot(third_trace)
plt.show()

与以前的模型不同，在这种情况下，我们的参数的后验分布是单峰的。

让我们总结一下这些后验分布。

pm.summary(third_trace)

除了查看汇总的统计数据，我们还可以查看地图解决方案。

**with** third_model:
    map_solution=pm.find_MAP()
d=dict()
**for** item **in** map_solution.keys():
    d[item]=[float(map_solution[item])]

third_map_coeffs=pd.DataFrame.from_dict(d)    
third_map_coeffs

我们看到，第二种和第三种模型的地图解决方案之间存在显著差异。预测呢？

**with** third_model:
    ppc = pm.sample_posterior_predictive(third_trace, samples=15000)*#compute y_score* 
**with** third_model:
    third_y_score = np.mean(ppc['y'], axis=0)
*#convert y_score into binary decisions* 
third_model_prediction=[1 **if** x >0.5 **else** 0 **for** x **in** third_y_score]*#compute confussion matrix* 
third_model_confussion_matrix =confusion_matrix(df['Occupancy'], third_model_prediction)
third_model_confussion_matrixarray([[1639,   54],
       [   3,  969]])

这个混淆矩阵与前两个模型中的相同。第二个和第三个模型的 y 分数分布如何？

figsize(12.5,4)
plt.hist(third_y_score-second_y_score, bins=40)
plt.title('Distribution of the difference between y_score and second_y_score')
plt.ylabel('Frequency')
plt.show(

这种分布集中在 0 左右。换句话说，y_scores 在不同模型间的分布几乎相同。系数呢，比如温度系数？

figsize(12.5,4)
plt.hist(third_trace['Temperature'][-40000:]-second_burned_trace['beta_temp'], bins=40, density=**True**)
plt.title('Difference between the temperature coefficients for the second and the third model')
plt.show()

差值近似遵循具有小平均值的正态分布。

让我们也检查一下湿度系数之间的差异。

figsize(12.5,4)
plt.boxplot(third_trace['Humidity'][-40000:]-second_burned_trace['beta_humid'])
plt.title('Difference between the humidity coefficients for the second and the third model')
plt.show()

同样，这种差异很小。

3.结论

我们看到，即使我们的模型使用不同的先验，预测性能是相似的。这证实了我们的信念，随着我们的数据集变大，它们应该收敛到同一个解决方案。

我们希望我们的项目将帮助 PyMC3 的初学者学习它的语法。我们发现 PyMC3 的代码相当直观，我们希望我们的代码清楚地证明了这一点。

4.参考文献。

[1]https://docs.pymc.io/notebooks/GLM-logistic.html官方 PyMC3 文档

[2]https://Goldin locks . github . io/Bayesian-logistic-regression-with-pymc 3/

PyTorch 上的贝叶斯 LSTM—使用 BLiTZ，PyTorch 贝叶斯深度学习库

原文：https://towardsdatascience.com/bayesian-lstm-on-pytorch-with-blitz-a-pytorch-bayesian-deep-learning-library-5e1fec432ad3?source=collection_archive---------22-----------------------

是时候给你的时间序列预测画一个置信区间了——现在这很简单。

LSTM 细胞插图。来源访问于 2020 年 4 月 14 日

这是一篇关于如何使用 BLiTZ 的帖子，BLiTZ 是一个 PyTorch 贝叶斯深度学习库，使用其贝叶斯 LSTMs 的实现来创建、训练和执行序列数据上的变分推理。

你可以查看这篇文章的示例部分的笔记本这里和 PyTorch 上 BLiTZ Bayesian 深度学习的知识库这里。

为了做到这一点，我们将解释贝叶斯长短期记忆是如何工作的，然后使用 Kaggle 的数据集通过一个股票置信区间预测的例子。

如果你是贝叶斯深度学习主题的新手，你可能想在 Medium 上寻找关于它的许多帖子中的一个，或者只是我们 lib repo 中关于贝叶斯 DL 的文档部分。你可能也想看看这篇关于 BLiTZ 用法教程的文章。

贝叶斯 LSTM 层

正如我们所知，LSTM 架构旨在解决当标准递归神经网络用于处理长序列数据时出现的信息消失问题。

在数学上，我们将 LSTM 架构翻译为:

LSTM 方程。来源:https://en.wikipedia.org/wiki/Long_short-term_memory由我用乳胶改写。访问时间为 2020 年 4 月 14 日。

我们还知道，贝叶斯神经网络的核心思想是，我们可以对它们进行采样以获得概率分布，然后优化这些分布参数，而不是具有确定性的权重。

利用这一点，有可能衡量预测的可信度和不确定性，这与预测本身一样，都是非常有用的洞察数据。

从数学上来说，我们只需要在上面的方程中增加一些额外的步骤。它们是权重和偏差采样，发生在前馈操作之前。

代表在层/模型的位置 N 上在第 I 时间采样的权重的方程。

代表在层/模型的位置 N 上第 I 次采样的偏差的方程。

当然，我们的可训练参数是参数化每个权重分布的ρ和μ。BLiTZ 有一个内置的 BayesianLSTM 层，可以为您完成所有这些艰苦的工作，因此您只需担心您的网络架构和培训/测试循环。

让我们来看看我们的例子。

首先，我们的进口

除了我们常见的进口产品，我们还将从blitz.modules引进BayesianLSTM，从blitz.utils引进variational_estimator一名室内装潢师，对他们进行不同的培训和复杂的成本收集。

我们还引入了collections.deque用于时间序列数据的预处理。

数据预处理

我们现在将创建并预处理数据集，以将其提供给网络。我们将从从 Kaggle 获得的数据集导入亚马逊股票定价，获得它的“收盘价”列，并将其规范化。

我们的数据集将由标准化股票价格的时间戳组成，并将具有形状(batch_size、sequence_length、observation_length)。

让我们导入并预处理数据:

我们还必须创建一个函数来转换时间戳中的股票价格历史。为此，我们将使用一个最大长度等于我们正在使用的时间戳大小的 deque。我们将每个数据点添加到 deque 中，然后将其副本附加到主时间戳列表中:

创建我们的网络类

我们的网络类接收到了variational_estimator装饰器，它简化了对贝叶斯神经网络损失的采样。它将有一个贝叶斯 LSTM 图层，其中 in_features=1，out_features=10，后跟一个 nn。Linear(10，1)，输出股票的标准化价格。

正如你所看到的，这个网络是一个非常正常的网络，这里唯一不常见的是BayesianLSTM层实例和variational_estimator装饰，但是它的行为是一个正常的火炬。

完成后，我们可以创建我们的神经网络对象，分割数据集，并继续进行训练循环:

创建对象

我们现在可以创建我们的损失对象，神经网络，优化器和数据加载器。请注意，我们没有随机分割数据集，因为我们将使用最后一批时间戳来评估模型。因为我们的数据集就大小而言非常小，所以我们不会为训练集创建数据加载器。

我们将使用正态均方误差损失和学习率=0.001 的 Adam 优化器

火车环线

对于我们的训练循环，我们将使用variational_estimator添加到我们的神经网络的sample_elbo方法。它对 X 个样本的损失进行平均，帮助我们轻松地进行蒙特卡罗估计。

为使该方法有效，网络的forward方法的输出必须与将被输入损失对象/标准的标签具有相同的形状。

评估模型并收集置信区间

我们将首先用要绘制的真实数据创建一个数据框架:

为了预测置信区间，我们必须创建一个函数来预测相同数据的 X 次，然后收集其平均值和标准差。同时，在参考真实数据之前，我们必须设置我们将尝试预测的窗口大小。

让我们看看预测函数的代码:

和置信区间的收集。请注意，我们可以决定距离均值有多少标准差来设置我们的置信区间:

由于我们使用的样本数量非常少，因此我们用高标准偏差对其进行了补偿。我们的网络将尝试预测 7 天，然后查询数据:

这里我们可以通过查看真实值是低于上界还是高于下界来检查置信区间。设置好参数后，您应该有一个 95%左右的置信区间，如我们所知:

查看预测图表

我们现在只是绘制预测图，直观地看看我们的训练是否顺利。我们将绘制真实数据和测试预测及其置信区间:

基于置信区间的股票价格预测。

为了结束我们的评估，我们将放大预测区域:

带有置信区间的测试数据的网络预测图。

结论

我们看到 BLiTZ 贝叶斯 LSTM 实现使得利用贝叶斯深度学习的所有能力在时间序列上实现和迭代变得非常容易。我们还看到贝叶斯 LSTM 很好地集成到 Torch 中，并且易于在任何工作或研究中使用和介绍。

我们还可以非常准确地预测 IBM 股票价格的置信区间，这可能是一个比点估计有用得多的洞察力。

参考

[## piEsposito/blitz-贝叶斯深度学习

BLiTZ 是一个简单且可扩展的库，用于创建贝叶斯神经网络层(基于 Weight…

github.com](https://github.com/piEsposito/blitz-bayesian-deep-learning) [## 了解 LSTM 网络

2015 年 8 月 27 日发布人类不是每秒钟都从零开始思考。当你读这篇文章时，你…

colah.github.io](https://colah.github.io/posts/2015-08-Understanding-LSTMs/) [## BLiTZ——py torch 的贝叶斯神经网络库

火炬动物园中的贝叶斯层是一个简单且可扩展的库，用于在…

towardsdatascience.com](/blitz-a-bayesian-neural-network-library-for-pytorch-82f9998916c7) [## 神经网络中的权重不确定性

我们介绍了一个新的，有效的，原则性和反向传播兼容的算法学习概率…

arxiv.org](https://arxiv.org/abs/1505.05424)

贝叶斯元学习

原文：https://towardsdatascience.com/bayesian-meta-learning-fa4eedcb89fc?source=collection_archive---------23-----------------------

这个故事介绍了贝叶斯元学习方法，其中包括贝叶斯黑盒元学习，基于贝叶斯优化的元学习，MAMLs 集成和概率 MAML。这是课程的简短总结‘斯坦福 CS330:多任务和元学习，2019 |第 5 讲——贝叶斯元学习’。

via 艺术字

这个故事涵盖了我们为什么需要贝叶斯方法的原因，以及如何实现和评估这些方法。它是课程‘斯坦福 CS330:多任务和元学习，2019 |第五讲——贝叶斯元学习’的总结。

对于元学习算法，3 个算法特性是重要的:表达能力、一致性和不确定性意识。 表现力是 f 代表一系列学习过程的能力，它衡量可扩展性和对一系列领域的适用性。一致性意味着学习的学习过程将解决具有足够数据的任务，这一特性减少了对元训练任务的依赖，从而导致良好的分布外性能。不确定性意识是在学习过程中对模糊性进行推理的能力。它让我们思考如何在强化学习的环境中探索新的环境，以减少我们的不确定性。它还考虑了如果我们处于安全关键环境中，我们想要校准不确定性估计。它也允许我们从元学习的贝叶斯观点来思考，什么样的原则方法可以从那些图形模型中导出？

这个故事涵盖了 1。为什么是贝叶斯？2.贝叶斯元学习方法 3。如何评价巴伊亚人

为什么是贝叶斯？

多任务原则和元学习原则是训练和测试必须匹配，任务必须共享“结构”。“结构”是什么意思？这意味着对共享潜在信息θ的统计依赖性。

来源:课程幻灯片

图形模型，其中θ是所有任务共享的肉参数。φI 是每个任务 I 的特定于任务的参数。特别是，这个图形模型显示了对θ的一些依赖性。基本上，每个参数ϕi 都有一个从θ开始的箭头。如果你的条件以θ为信息，任务参数(ϕs)就变得独立，也就是ϕi 独立于另一个以θ为条件的ϕi。如果你把θ作为潜在信息的条件，那么它们就不是独立的。这些性质表明给定θ的ϕ上的分布比ϕ.的边缘分布具有更低的熵本质上，θ告诉你关于特定任务参数ϕ.的信息

如果您可以识别θ(例如，通过元学习)，**在哪些情况下，学习ϕ不会比从头开始学习更快？如果ϕs 在知道以θ为条件之前是独立的，这意味着它们的熵将是相同的，那么从零开始学习将与从共享信息θ中学习一样快。如果单个点承载了θ的所有信息，那么共享的信息基本上不会告诉你那么多，因为信息也存在于单个数据点中。以θ为条件的熵会相当高，边际熵也会相当高。**如果给定θ的ϕ熵为零会怎样？这意味着在ϕ中没有任何附加信息不被θ捕获，那么θ可以解决所有的任务，你不需要做任何事情来学习θ。

【theta 可能包含什么信息？如果任务族对应于具有不同振幅和不同相位的正弦曲线，那么θ将对应于正弦函数族。在机器翻译的例子中，θ对应于所有语言对的家族，而ϕ中没有出现在θ中的信息将对应于特定语言的特定事物。注意，在这两个例子中，θ比所有可能函数的空间都要窄(这就是为什么原则上我们可以从使用元学习中获益)

如果你有无限多的任务，你应该可以精确地恢复θ，或者基本上高精度地恢复那一族。如果你在没有大量任务的情况下进行元学习会怎样？如果你有一个任务空间，它不一定会覆盖任务的真实分布，但可能会对任务空间进行元过拟合，这样实际上不会恢复对应于所有语言对的θ，但它会找到对应于一组语言对的θ，这些语言对看起来像你的训练数据中的东西，而不是捕捉完整分布的东西。因此，除非测试数据非常接近训练样本，否则你不会有效地适应来自该分布的新事物。

少数几个学习问题可能是模糊的。有“+”和“-”两类。对于“+”类，图片的特征可能是“微笑”、“戴帽子”和“年轻”。但是最右边的测试图像包含一些特征，但不是全部，所以很难区分这些测试图像属于哪一类。来源:课程幻灯片

参数方法使用这种分布的确定性估计，而贝叶斯元学习方法试图从这种分布中采样。来源:课程幻灯片

对于**参数方法，**在给定数据集和元参数θ的情况下，我们恢复的是任务特定参数ϕ的确定性估计。所以你会得到这个分布的点估计。这在哪里/为什么会成为问题？在某些情况下，我们需要的不仅仅是点估计。例如，一些少量的学习问题可能并不完全由他们的数据决定，也许在你有先验信息的证据下，基础函数是不明确的。在本例中，模糊性是因为数据集很小。为了解决这个问题，**如果我们可以生成关于底层函数的假设，会怎么样？**如果我们可以从这个分布中取样，我们就可以对我们的不确定性进行推理，这对于安全关键的少量拍摄学习非常重要。它允许我们在元学习中主动学习和探索。这就是贝叶斯元学习方法的用武之地。

贝叶斯元学习方法

从计算图的角度来看。来源:课程幻灯片

如果我们关心在我们的预测上生成分布， 一种方法是让函数 f 输出在标签上分布的参数 。例如，如果你的标签 y 是离散的，你让 f 输出对应于离散分类分布概率的东西；您还可以输出高斯分布的均值和方差来表示标签上的高斯分布；或者，如果你有一个更多模态的分布，你可以试着表示高斯混合的均值、方差和混合权重；或者，如果您有一个多维输出标签 y，那么您可以输出一系列条件分布的参数，以允许您表示这些变量的联合分布，这是一个自回归模型。一旦你输出了分布的参数，你就可以用 最大似然估计优化 来训练这些方法，优化分布的均值和方差，或者基本上优化该分布的所有全概率值。

这种方法的好处是它们简单并且可以与多种方法结合。缺点是这种方法不能推理底层函数的不确定性，以确定数据点之间的不确定性如何相关。它还倾向于**产生校准不良的不确定性估计。**例如，你不能把模型不确定性和标签噪声分开。它也被限制在 y 上的一类特殊分布。

我们能为ϕ做同样的最大似然训练吗？我们不能以完全相同的方式来做，因为我们没有对应于ϕ的标签，除非我们生成标签。

贝叶斯元学习方法

贝叶斯深度学习工具箱

一种可变自动编码器。来源:课程幻灯片

贝叶斯深度学习旨在用神经网络来表示分布。有许多方法可以用神经网络来表示分布。 一种流行的方法 就是使用潜变量模型，然后用变分推理对其进行优化。 具体来说，我们有一个图形模型，其中有潜在变量 Z 和观察变量 x。例如，左边的图形模型是一个可变自动编码器。你可以制定一个似然目标的下限，并使用它来优化 X 上的分布，其中 X 上的分布可能是非高斯的，因为它有这个潜在变量。

另一种方法在某些方面非常流行，也非常简单，那就是 使用基于粒子的分布图 。特别是，您可以在不同的数据引导上训练不同的模型。这些模型中的每一个都对应于你的一个特定的分布，然后你把这些粒子组合起来，代表那个分布的样本。

另一种方法 是表示神经网络参数的权重的显式分布。 然后练习这些分布趋向于带有对角协方差矩阵的高斯分布，这样你基本上就有了每个神经网络参数的独立方差。这使您可以表示函数的分布，但需要注意的是，这种两个参数独立的独立假设在实践中总是被违反。

规格化流程 通过将一些潜在分布转化为你的数据分布，或者从潜在空间转化为你的数据空间，再转化回你的潜在空间，尝试表示数据分布上的函数 。

在 e 基于能量的模型和 GAN 、 中，估算非标准化密度的方法是对数据使用低能量，对其他所有数据使用较高能量，而其他所有数据都由 GAN 中的生成器近似。

埃尔博。来源:课程幻灯片

这个图形模型包含了观察变量 x，和潜在变量 z .T3【证据下界】或者 ELBO 希望能够估计出一个关于数据的似然性的下界 T5【以便优化数据的似然性。ELBO 可以由关于 q(z|x)的期望加上在 q(z|x)上正则化或操作的熵项来表示。q 是表示 x 和 z 的概率的变分分布，这个可以改写为两项互相加， 第一项对应的是解码器的重构损失，基本上是从推理网络 q 采样后根据解码器的数据的似然性，第二项对应的是推理网络和先验之间的 KL 散度。 P 对应于模型，q 对应于为了逼近似然目标 而引入的变分分布。对于如何表示 p，有几种设计选择:p(z)和 p(x|z)。 P (x|z)可以表示为一个神经网络， 例如，在变分自动编码器的情况下，它表示一个神经网络，该网络接收潜在变量并输出图像或您正在建模的任何数据类型。P(z)通常被表示为具有单位方差的对角高斯。 它也可以由神经网络来表示，或者表示为学习的平均值和学习的方差。实际上，变分自动编码器通常不学习它，因为之后的层可以将它转换成学习的均值和方差。 q(z|x)也是用神经网络来表示的，推断为推理网络或者你的变分分布。 在深度学习的上下文中使用变分推理时，θ表示模型参数，ϕ表示推理网络的参数。

重新参数化的把戏。来源:课程幻灯片

但是优化 ELBO 有个问题。ELBO 包含关于 q 的期望，这意味着我们需要能够反向传播到 q 分布中，但是采样是不可微分的 e。我们可以使用 重新参数化技巧 。*对于高斯 q(z|x ),它可以用一种可微分的方式来表示，*用从单位高斯采样的噪声来重新参数化。特别是，如果我们将潜在变量的分布表示为与输出均值和方差的神经网络 q 的输出相对应的高斯分布，那么您可以将 q 的输出(高斯分布)表示为由噪声重新参数化，而不是从均值加上方差乘以噪声的特定分布中采样。幸运的是，这个方程是可微的。如果你的推理网络足够有表现力，它应该能够将你的数据分布转换成潜在变量的高斯分布。

这通常被称为 摊销变分推断 ，我们有一个推断网络，它预测变分分布，基本上摊销了估计该分布的过程。元学习可以用摊销变分推理吗？

贝叶斯元学习方法

贝叶斯黑盒元学习

黑盒元学习方法。来源:课程幻灯片

首先，为了简单起见，让我们考虑一下黑盒元学习方法。特别是，我们希望有一个神经网络，它将训练数据集作为输入，并在我们的参数φ上产生分布。然后应用ϕ对神经网络进行参数化，该神经网络将 x 作为输入，将 y 作为输出。目标是在根据我们的训练数据集对特定任务参数进行采样时，能够最大化给定特定任务参数的测试数据点的可能性。

摊提变分推理成元学习。来源: 课程幻灯片

在标准 VAE 中，观察变量是 x，潜在变量是 z。而在元学习中，观察变量是数据集，潜在变量是ϕ.ELBO 开始学习数学。它在特定任务参数ϕ上有一个变化的分布，在这里我们将从给定的ϕ.中取样并估计数据的可能性 KL 项表明 phi 上的变分分布和先验应该是相似的。

来源:课程幻灯片

对于 q 条件有两种设计选择。如果您希望能够将参数作为我们数据集的函数进行采样，那么我们应该对我们的训练数据进行条件 q(这显示在第一个等式中)。这个训练数据和测试数据点的概率有什么不同？ 因此，当根据我们的训练数据集 对我们的任务特定参数进行采样时，您希望能够最大化给定我们的任务特定参数的测试数据点的可能性(如第二个等式所示)。

元参数θ。来源:课程幻灯片

元参数θ呢？我们将有一个以θ为条件的特定任务参数ϕ的先验。测试数据点的分布仅是φ的函数。我们也可以θ为条件。

最终目标(为了完整性)。来源:课程幻灯片

最终目标是使所有任务的元参数和期望值最大化。特定于任务的参数从神经网络 q 中采样，然后应用它们来最大化 p 的可能性。右边的一项鼓励 q 接近某个先验分布。

贝叶斯黑盒元学习概述。来源:课程幻灯片

为了总结贝叶斯黑盒元学习，列出了一些好处和缺点。一个好处是它可以表示 y 上的非高斯分布，并且它还**生成任务特定参数ϕ上的分布，**而不仅仅是生成标签上的分布(y^{ts}).这使得它能够表示基础函数的不确定性，而不仅仅是基础数据点的不确定性。缺点是在给定θ的情况下，它只能表示ϕp 的高斯分布。这有两个原因。一个是重参数化技巧，它适用于高斯。第二件事是，在目标中出现的 KL 项，对于高斯型目标是可以以封闭形式评估的，但是对于其他非高斯型目标是不能以封闭形式评估的。

贝叶斯元学习方法

基于贝叶斯优化的元学习

将基于梯度的元学习重铸为分层贝叶斯。来源:课程幻灯片

我们可以将基于梯度的元学习重铸为分层贝叶斯。这很好，因为它将 MAML 的贝叶斯解释作为一种学习元参数，这样在测试时，您可以在元参数所代表的高斯先验下进行地图推断。它使用特定于任务的参数ϕ的地图估计作为θ的函数，这意味着它表示该分布的点估计，并且它只给出该分布的一组参数，例如，您无法从该分布中进行采样或更容易地表示任务参数的完整分布。

贝叶斯元学习方法

摊销贝叶斯元学习

Q 是任意函数 ，它可以包含一个梯度算子。这是*，所以你可以有一个推理网络，基本上是在你的推理网络*内执行梯度下降。因为 q 不必是神经网络，它可以是输出ϕ.分布的任何东西 (摊销贝叶斯元学习) 在文中，我们将运行梯度下降 w.r.t .一组参数的均值和该组参数的方差 w.r.t .一些训练数据，然后其他一切都是一样的。 与 MAML 相比，这不仅仅是对参数进行梯度下降，*这是对参数的均值和方差进行梯度下降，在这个过程的最后，你会得到参数的均值和方差，而不仅仅是均值。好处是你在测试时运行*梯度下降，**所以我们得到一个基于优化的元学习方法，基本上是把梯度下降填充到参考网络中。缺点是 p(ϕi|θ)被建模为高斯。

我们可以建立非高斯后验模型吗？

贝叶斯元学习方法

合奏

我们可以通过维护单独的模型实例来模拟非高斯后验。有几篇关于组装哺乳动物的论文。

MAMLs(aka EMAML)的集成训练 M 个独立的 MAML 模型，如果集成成员太相似，就不能很好地工作。我们可以将这种想法应用到集成黑盒和非参数方法中。

来源:课程幻灯片

BMAML 是一个更加多样化的哺乳动物群体，基本上，它将不同的群体成员相互推开，鼓励他们代表不同的模型。 实现这一点的方法是，当你运行梯度下降时，你有典型的似然项，但你也有一个鼓励不同系综成员彼此不同的项。你可以用不同的核函数来表示模型之间的相似性。

BMAML 优化。来源:课程幻灯片

他们不仅仅是将粒子相互推开，还进行了优化，使集合给你一个产生高可能性的粒子分布。因此，他们没有完全分别训练这些不同的模型，而是把它们组合在一起，优化了一个术语，共同优化了它们。所以你在优化这 M 个粒子的平均可能性。

这些集成模型的好处是它们简单，往往工作良好，并且它们可以模拟非高斯后验。缺点是您需要维护 M 个模型实例，或者只在最后一层进行基于梯度的推理。我们可以在所有参数上模拟非高斯后验模型而不必维护单独的模型实例吗？

贝叶斯元学习方法

概率 MAML

我们可以用像哈密尔顿蒙特卡罗这样的程序对参数向量进行采样吗？具体来说，哈密顿蒙特卡罗所做的是添加噪声，然后重复运行梯度下降，以便能够从一些分布中采样。直觉是学习一个先验，一个随意的踢腿可以让我们进入不同的模式。

损失景观。来源:课程幻灯片

在这种损失格局中，有不同的解决模式。一个对应于分类“微笑”和“戴帽子”的分类器，一个对应于分类“微笑”和“年轻”的分类器。我们希望参数向量在中间，这样如果我们给参数向量添加一些噪声，然后运行梯度下降，我们得到这个解的两种不同模式。我们得到不同的函数，代表正确答案的不同模式。

概率 MAML 中测试时的推断。来源:课程幻灯片

在这种情况下，θ将是一个分布，而不是一个单一的参数向量。我们的目标是在测试时给定所有观察到的变量，即x train、y train和x test，对不同的特定于任务的参数进行采样。我们可以跨越x test，因为y test没有给定，其他的保持条件独立，所以你可以从ϕ给定的函数x train和y train中采样。这种分布的精确解是完全难以处理的，因为它需要积分。如果我们知道如何在给定θx train和y train的情况下对 phi 进行采样，那么采样就会变得容易得多。通过地图推断，我们可以用ϕ点估计来粗略地近似这个分布。实际上，我们仍然可以从ϕ采样，首先从θ采样，然后我们运行几个梯度下降步骤。这是在测试时我们想做推理时发生的。训练可以用摊余变分推理来完成，详情可以查阅论文。

概率 MAML 中的祖先抽样。来源:课程幻灯片

祖先取样是什么样子的？有 2 个步骤: 首先我们从θ的 P 中采样，这个 P 有均值和方差；然后我们运行梯度下降，从样本θ开始，运行梯度下降到θ。它允许我们表现这些类型的多模态分布。

这种方法的好处是它具有非高斯后验，在测试时简单，并且只有一个模型实例(与 MAML 的系综相反)。缺点是它有更复杂的训练程序。

贝叶斯元学习方法

方法概述

涵盖了三种方法。第一个是训练函数 f 来输出在y test上的分布。好处是它很简单，可以结合多种方法。但是它不能推理基础函数上的不确定性，并且也可以表达y test上的有限类分布。第二种方法是黑盒方法，其中使用潜在变量模型和摊销变分推断。好处是它可以代表y test上的非高斯分布，但它只能代表高斯分布 p(ϕi |θ)，当ϕi 是一个潜在向量时是可以的。对于基于优化的方法，我们提到了 3 种类型:摊提推理；集成和混合推理。摊销推断很简单，但它将 p(ϕi |θ)建模为高斯分布。集成方法简单，往往工作良好，具有非高斯分布，但它包含 M 个模型实例。混合推理具有非高斯后验，测试时简单，只有一个模型实例，但训练过程比较复杂。

如何评价巴伊亚人

如何评价一个贝叶斯元学习者？

使用标准基准怎么样，比如 Minilmagenet 精度？好处是它是标准化的，所以我们可以更容易地比较论文。Minilmagene 也有实像。但是像准确性这样的度量标准不能评估不确定性。任务可能不会表现出模糊性，并且不确定性在该数据集上可能没有用。什么是更好的问题和度量？这取决于你关心的问题。

检查这些算法性能的一种方法是看不同的玩具例子。例如，设想您从中取样的基础函数。

来源:课程幻灯片

这个例子对应于基本的元学习函数是正弦函数和线性函数，其中标签中有噪声，使得任务不明确。它所推理的函数类型，在某些情况下，是模糊的，即使它是正弦函数或者线性函数。你可以从抽样分布中看到这种多模态输出。这表明函数实际上可以表示多模态的东西。

来源:课程幻灯片

第二个例子是一个分类的例子，你只给了一个数据点，所有的任务都对应于循环决策边界。你试图区分积极点和消极点。D train仅包含一个正数据点，而D test 包含正负数据点。因此，基本上你可以看到虚线所示的决策边界的可视化，其中神经网络代表这些不同的决策边界，这些边界代表了它接受训练的函数类的结构。

我们也可以使用可能性或准确性的组合以及分布的覆盖范围来衡量性能。或者我们可以绘制预测器的置信度与预测器的精确度的关系图。还有主动学习评价。

参考

1. 斯坦福 CS330:多任务和元学习，2019 |第 5 讲——贝叶斯元学习
2. 课程幻灯片
3. 具有潜在嵌入优化的元学习
4. 元数据集:用于从少量示例中学习的数据集的数据集
5. 主动一次性学习—斯坦福计算机科学
6. 主动学习的学习算法
7. 自动编码变分贝叶斯
8. 使用深度集成进行简单和可扩展的预测不确定性估计
9. 概率模型不可知元学习
10. 概率模型不可知元学习——NIPS 会议录

英国停车搜索的贝叶斯模型

原文：https://towardsdatascience.com/bayesian-modelling-of-uk-stop-and-search-3f9947272fce?source=collection_archive---------28-----------------------

用贝叶斯推理拟合泊松广义线性混合模型研究停止和搜索的过度使用

国际国王教会在 Unsplash 上拍摄的照片

在英国和国际上，关于过度使用拦截和搜查的讨论越来越多。看看英格兰和威尔士警察部队拦截和搜查次数与逮捕次数的比例，我们可以调查不同种族群体之间是否存在差异。具体来说，泊松广义线性混合模型被实施以允许种族群体固定效应和警察力量随机效应。该模型符合使用 Gibbs 抽样的贝叶斯推断，并且我们发现在样本中对于我们的种族群体效应的实质性差异。这证明了警察部队在地理层面过度使用拦截和搜查。这个项目的 bug 模型的实现可以在 GitHub 上找到。

背景

现有文献旨在调查过度使用拦截和搜查。Bowling 和 Phillips (2007 年)提供了关于英国文献的广泛总结和讨论，认为统计数据显示对黑人个人使用拦截和搜查权力是不成比例的。这里要看的一个常见统计数据是按种族划分的每个居民的停留和搜索次数。这可以从折线图中看出，例如，在 2019 年，每 1000 名黑人居民中黑人个体被拦下和搜查的次数比每 1000 名白人居民中白人个体被拦下和搜查的次数多 9 倍以上。

Bowling 和 Phillips (2007)强调，按常住人口标准化的统计数据经常受到批评，因为没有考虑潜在的不同犯罪参与和“可用性”。在这里，Miller (2000)将“可用”定义为描述一个人在停车和搜查发生的时间和地点使用公共场所。此外，Miller (2000)确实发现了常住人口和可居住人口的差异，与常住人口相比，来自黑人、亚洲人和少数民族背景的个人在可居住人口中的比例过高。因此，在这个项目中，拦截和搜查是按逮捕比例进行调查的。但是，应该注意的是，如果在拦截和搜查过程中存在偏见，在逮捕过程中也可能存在偏见，因为这两个过程通常都需要“合理的理由”。似乎可以合理地假设，如果存在的话，这些偏见将具有相同的方向，因此，按照逮捕的比例来看待拦截和搜查可以作为不相称的最低门槛。

数据

拦截和搜查个人层面的数据是公开的，由英国警方提供。对于这个项目，我们已经采取了 2015 年 1 月至 2020 年 2 月(含)英格兰和威尔士警方报告的所有拦截和搜查(约 190 万次)。探索这些数据很有趣，例如，最常见的搜索对象是毒品(56%)，而攻击性武器仅占所有搜索对象的 10%。此外，大约 70%的拦截和搜查导致没有采取进一步行动，只有大约 15%导致逮捕。如上面的堆积条形图所示，除了逮捕之外，还可以采取范围广泛的进一步行动。

每份拦截和搜查报告都包括警察力量和自我定义的种族。因此，拦截和搜查次数是按种族和警察队伍汇总的。由 Office⁴.内政部提供的按族裔群体分列的被捕人数也是公开的这些是在 2015 年 3 月至 2019 年 3 月期间收集的。该地图显示了“拦截搜查差距”，即黑人被拦截和搜查的次数与黑人被逮捕的次数的比例比白人被拦截和搜查的次数与白人被逮捕的次数的比例。停止逮捕差距大于 1 反映了黑人相对于白人更频繁地被停止和搜查，与逮捕成比例。警察力量似乎有很大差异，支持在我们的泊松广义线性混合模型中使用随机警察力量效应。

泊松广义线性混合模型

泊松广义线性混合模型的实施考虑到了种族群体的固定效应和警察力量的随机效应，可以被看作是盖尔曼和 al.(2007)⁵.的扩展作者用种族群体固定效应和警察辖区随机效应研究了纽约的拦截和搜查。因此，该项目的实施着眼于更大范围的警察部队与辖区的区域单位以及英格兰和威尔士与纽约的地理区域。我们的最终模型如下所示。

该模型采用马尔可夫链蒙特卡罗方法进行贝叶斯推理拟合。具体来说，选择 Gibbs 抽样是因为模型的层次结构，并且使用了 BUGS 软件。假设所有未知参数都有非信息先验。运行三个平行链，并结合轨迹图使用 Gelman-Rubin 统计来评估收敛性。通过 2，000，000 次迭代和 1，000，000 次老化，达到了令人满意的收敛证据。10 的细化用于适应 bug 内存限制。

结果

结果表显示了通过取样本平均值估计的后验平均值。此外，可信区间代表 95%的高后验密度区间。任何种族群体截距的可信区间都不会跨越零，因此给定种族群体效应的值大于零的概率为 95%。由于对数关联函数，必须采用后验均值的指数来量化族裔群体对警察部队拦截和搜查次数的影响，以及与逮捕的比例。因此，通过取后验平均值的指数，我们可以看到，亚裔、黑人和混血儿被拦截和搜查的次数分别是白人的 1.52 倍、1.30 倍和 1.20 倍，与被逮捕的人数成比例。最后，下面的轨迹图显示在左侧，支持收敛的证据。此外，在右侧，种族群体效应的密度图看起来平滑且呈单峰。此外，零度以下的密度非常低。

结论

由于难以控制潜在变化的犯罪参与和可利用性，衡量过度使用拦截和搜查存在很大困难。该项目通过观察拦截和搜查与逮捕的比例来解释这些情况，作为不相称的最低门槛。将警察部队一级的拦截搜查和逮捕汇总，发现不同族裔群体的拦截搜查率与逮捕率之间存在差异。据估计，亚洲人、黑人和混血儿被拦截和搜查的频率分别是白人的 1.52 倍、1.30 倍和 1.20 倍，与被捕人数成比例。这与 Gelman 等人(2007)在纽约警区的发现一致。有缺陷的 Gibbs 抽样被证明是模拟后验分布的一种灵活有效的方法。此外，随机效应的方差参数的可信区间不跨越零，支持混合效应模型的实现。最后，如果在较低的地理水平上可以获得逮捕数据，那么用较小的区域单位重复这种分析，以调查是否有证据表明警察部队内部的随机效应，将是有用的。

参考

[1]鲍林，b .和 c .菲利普斯(2007)。不相称和歧视:审查警察拦截和搜查的证据。现代法律评论。

[2]米勒，J. (2000 年)。可用于停靠和搜索的人口分析。内政部警察研究系列。

[3]英国警方(2020 年)。停下来搜索数据。【https://data.police.uk/data/】T4。

[4]家庭办公室(2019)。警察权力和程序英格兰和威尔士统计。https://www . gov . uk/government/collections/police-powers-and-procedures-England-and-Wales。

[5]赫尔曼、法甘和基斯(2007 年)。在种族偏见的背景下对纽约市警察局拦截搜身警察的分析。美国统计协会杂志。

贝叶斯网络和合成节点

原文：https://towardsdatascience.com/bayesian-networks-and-synthetic-nodes-721de16c47e2?source=collection_archive---------37-----------------------

之前，我已经直观地介绍了什么是 bn以及如何构建它们。在这篇文章中，我给出了一个更正式的观点，并讨论了建模的一个重要元素:合成节点。

虽然神经网络基于可用数据给出答案，但贝叶斯网络(BNs)包括与我们的研究相关的非样本或先验人类专业知识。通常，采访专家询问几个参数的影响是收集所有信息的最佳方式，这使得贝叶斯网络比其他技术更丰富。将这种专业知识与样本数据相结合使这种技术变得强大。此外，他们提出了模型的图形表示结构，在这种情况下，这比神经网络等黑盒模型更好。更容易解释获得的结果并与专家交流:当与不是建模技术专家而是商业领域(即医疗保健或风险管理)专家的人一起工作时，这提供了优势。BNs 的另一个优势是它们能够通过专家知识和原始数据来构建，即使先前的知识是不完整的。可以从专家的知识开始，用数据提炼，然后用专家的知识进行反馈，等等。

贝叶斯网络综述

贝叶斯网络(BN)是概率模型的一部分，该概率模型通过图表示一组变量及其依赖关系，即一对 G = (V，e)，其中 v 是不同节点的有限集合，E ⊆ V×V 是弧的集合。有序对( u，v ) ∈ E 表示从节点 u 到节点 v 的有向弧， u 被称为 v 的父节点， v 是 u 的子节点(Pearl 1988)。

BN 通常是表示概率因果关系的有向无环图(DAG ),其中节点表示变量，弧线表示变量之间的依赖关系。更准确地说，节点是随机变量的状态，从父节点指向子节点的弧线是两个节点之间的因果条件依赖。假设父节点状态的概率不同，则因果关系由节点状态的概率表示。

BN 使用条件概率来描述事件。任何关于事件或假设的不确定性的信念都是暂时的。这被称为先验概率或“P(H)”。这个先验概率被新的经验“e”更新，提供了关于“H”的不确定性的修正信念。给定证据的新概率称为后验概率或“P(H|e)”。贝叶斯定理描述了后验概率和先验概率之间的关系“P(e|H)”是假设假设为真，证据为真的概率，而“P(e)”表示在所有可能的假设下，新证据为真的概率。

条件概率表(CPT)展示了节点之间的概率关系。下图在三个基本图表上显示了这种关系。在这个例子中，每个节点有两种状态:真(T)和假(F)。该图还显示了每个节点配置的不同 CPT。

三个基本的 BN 模型

对于 BN，条件概率分布可以表示为一组参数:

其中 i = 1，…， n 定义每个节点； k = 1，…，r i 定义每个 r i 变量值(节点的状态)；并且 j= 1，…，qi 定义了父变量节点的 qi 变量配置的集合(例如上图中的每个图形配置)。因此，对于父节点的 q 个 i 个有效配置中的每一个，存在具有向量的多项式分布

例如，在上图的共同效果图中，黄色 P(Z=T | X=T，Y=F)中的参数表示假设 X 为真，Y 为假，Z 为真的条件概率。

BN 配置的参数化需要矢量的完全估计。随后，两种类型的分析是可能的:(1)反向推理，它允许给定的一个观察发现最可能的原因中的假设(诊断)，和(2)自上而下的推理，它允许估计概率的一个观察给定的假设(预测)(Pearl 1988)。

使用合成节点

在 BN 中，对于由 pa (Xi)表示的节点 Xi 的一组父节点，节点值的联合概率分布可以写成每个节点及其父节点的局部概率分布的乘积:

计算目标节点的概率分布所需的变量的数量随着父节点的数量及其状态的数量而增加，这是耗时的。合成节点应该减少计算目标节点 CPT 所需的数据量。正如(康斯坦蒂努等人，2016)提到的，“合成节点是一个简单地由它的父节点的值使用一些专家驱动的组合规则定义的节点”。因此，合成节点对于减少因果复杂性和 CPT 中组合爆炸的有害影响是有用的。此外，这些节点改善了 BNs 的结构，即因果关系的可视化。通过使用合成节点，CPT 的大小和复杂性“大大降低，提高了模型的计算速度和准确性(Fenton 和 Neil 2013)。

例如，上图(a 部分和 b 部分)显示了从简单网络开始使用合成节点的方案。在这个例子中，网络有四个变量 a、b、c 和 d，每个变量有四个状态，计算 P(A|D，b，c)的 CPT 的概率值的数目是 4^4 = 256。相反，通过引入合成节点 s，该计算可以简化为两个表 P(A|B，s)和 P(S|C，d)。现在，该模型需要使用 4^3 + 4^3 = 64 值而不是 256 值来计算 s 和 a 的两个 CPT。

(Neil 等人，2000)中描述的合成节点的使用示例

“从从业者的角度来看，考虑到当前算法的准确性和速度，使用最新的软件工具编译和执行 BN 的过程相对来说没有痛苦。然而，为特定问题构建完整的 BN 的问题仍然存在，即 如何构建图结构以及如何为图 的每个节点定义节点概率表。”(尼尔、芬顿和尼尔森 2000)。

感谢阅读！！

如果你想理解贝叶斯网络的逻辑(一个更基本的信息)，我推荐你这篇文章:

[## 贝叶斯思维导论:从贝叶斯定理到贝叶斯网络

假设世界上存在一种非常罕见的疾病。你患这种疾病的几率只有千分之一。你想要…

towardsdatascience.com](/will-you-become-a-zombie-if-a-99-accuracy-test-result-positive-3da371f5134)

如果你想了解构建贝叶斯网络的参数(和超参数)是什么，我建议你这样做:

[## 贝叶斯网络中的超参数调整问题

在这段历史中，我们讨论了在基于贝叶斯理论建立模型时要考虑的结构标准

towardsdatascience.com](/the-hyperparameter-tuning-problem-in-bayesian-networks-1371590f470)

参考书目

康斯坦蒂努、安东尼·科斯塔、诺曼·芬顿、威廉·马什和卢卡斯·拉德林斯基。2016."从复杂的问卷调查和访谈数据到用于医疗决策支持的智能贝叶斯网络模型."医学中的人工智能 67:75–93。

芬顿、诺曼和马丁·尼尔。2013.用贝叶斯网络进行风险评估和决策分析。Crc 出版社编辑。

尼尔、马丁、诺曼·芬顿和拉斯·尼尔森。2000."建立大规模贝叶斯网络."知识工程评论 15(3):257–284。https://doi.org/10.1017/S0269888900003039.

**珠儿，朱迪亚。1988."智能系统中的概率推理."摩根考夫曼圣马特奥。https://doi.org/10.2307/2026705.

贝叶斯神经网络:1 为什么烦恼？

原文：https://towardsdatascience.com/bayesian-neural-networks-1-why-bother-b585375b38ec?source=collection_archive---------33-----------------------

贝叶斯神经网络

由法鲁尔·阿兹米在 Unsplash 上拍摄

这是贝叶斯深度学习系列的第一章。下一章可用 此处 而第三章则是 此处 。

我经常被要求解释我对贝叶斯神经网络的古怪兴趣。所以这是一系列关于他们的开始。因此，我可以有把握地说，这是一系列关于深度学习未来的开始——一个现在就可以实现的未来。在接下来的几周里，您将学习如何使用为贝叶斯型神经网络构建的一个新库。我将包含代码并讨论 TensorFlow-Probability 和 Pytorches Pyro 的工作，包括安装、训练和预测各种不同的网络架构。在本章中，我们将讨论以下目标:

了解为什么贝叶斯神经网络如此有用和令人兴奋
了解它们与普通神经网络的不同
理解您可以获得的不确定性度量是一个主要优势
安装一个库，开始处理所有随机的事情

今天，许多深度学习的开发者并不是从计算机科学或学术背景开始的。反复试验，对一本书和一些博客的耐心大有帮助。不出所料，贝叶斯深度学习不需要任何不同。

但是一些基础理论知识的优势比平时更大。令人放心的是，这可以压缩，表达和理解相当简单。首先，是不寻常的术语引起了混淆。因此，在开始之前，我们将用一些代码来讨论它们。

我们对随机变量感兴趣。这意味着大多数变量在需要的时候都会改变值。所以非正式的随机变量是不可靠的东西(更正式的说法是我们称之为‘随机变量’，但是谁会从正式中受益呢？).使用张量流概率(TFP)或 Pytorches Pyro，您可以调用:

rv_normal = dist.Normal(loc=0., scale=3.)
rv_normal.sample([1])

这将返回一个值，每次我们运行第二行时，这个值都是不同的。即使我们不重新运行第一行也是如此，因为 rv_normal 是由正态(高斯)分布定义的随机变量。

因此，您可以将这些随机变量插入到神经网络中，并使用它们来代替我们通常用于权重的点值。你为什么要这么麻烦？因为在网络被训练之后，我们所做的每一个预测都将具有不同的随机可变权重。所以，整个模式会不一样。然而——在我放开你之前，我们用随机变量构建的神经网络的重要性在于，这些变量就像被拴在狗链上一样，只能返回特定范围内的值。这是因为我们的神经网络训练已经训练了我们的正态分布，使其钟形峰值位于最佳值附近。

从所有这些中可以得出三点:
分布拥有随机变量，每次我们做预测时，这些随机变量都会返回我们的随机权重值。
我们的神经网络训练移动整个分布，而不仅仅是一个点权重值。
这类似于在传统神经网络中使用分布初始化权重的过程(用于权重的初始值)。除了我们在整个训练和预测中继续使用分布，而不仅仅是初始化。

训练分布更正式地称为调节分布(类似于巴甫洛夫调节他的狗。(我似乎既不能逃避形式，也不能逃避狗)。所以所有这些复杂的随机变量和分布加起来就是一件事…，不同的预测。但是，如果我们已经很好地训练了贝叶斯神经网络，我们的预测将非常相似(对于连续、浮点、回归输出)和相同(对于分类输出)。).如果我们训练我们的神经网络不好，我们每次都会得到非常不同的结果。如果贝叶斯网络可以说话，它会说它的预测相似或非常不同之间的差异是它的不确定性。将不确定性视为信心(只有贝叶斯神经网络不喜欢谈论他们的信心，所以称它为不同的东西)。如果我们建立了一个良好的网络架构，最重要的是如果我们提供了良好的数据，神经网络将给出自信的预测，如果我们做的一切都不好，我们将得到不自信的预测。在贝叶斯魔法术语中，当我们给出坏消息时，神经网络会更加不确定，而当我们给出好消息时，神经网络会更加不确定。术语上的神奇差异因为贝叶斯神经网络从来都不是确定的，它们只是不太确定。而我们传统的神经网络总是确定一切！当你第一次开始使用神经网络时，你可能已经看到了类似的东西。您完成了对第一个模型的训练，并且表现一般，然后您在不做任何更改的情况下训练了另一个模型，平均来说仍然表现一般，只是您得到了完全不同的结果。贝叶斯神经网络就像那样，只是不是依靠你的挫折来保持训练有微小差异的神经网络，而是通过贝叶斯方法，你训练一次，就内置了无限数量的模型。

太棒了。掌声。那又怎样？所以…我们可以改变这一切，通过观察结果的可变性来看我们的数据(主要)有多适合这项任务。但比贝叶斯模型好得多的是，我们通常能够用比神经网络常规所需数据更少的数据来训练模型。贝叶斯模型将很好地推广。它也较少受到每个训练类的数据量不平衡的困扰——这是表现不佳的一个被忽视的原因。

说够了，我的好人，告诉我们怎么做。对哦。

首先，您需要安装两个主要概率编程库中的一个，以遵循本系列。虽然您可以仅使用核心库来构建贝叶斯神经网络，但您需要手动实现训练代码，这不会在数学和隐藏错误方面感到害羞。所以我们要么用 Pyro(建立在 Pytorch 之上)，要么用 TensorFlow-Probability(可能很明显建立在 TensorFlow 之上)。如果您有使用某个核心库的经验，最好继续使用它，并从您的熟悉中获益，而不是改变它。滚动到你喜欢的一个相关的库部分，或者如果你是一个受虐狂，继续执行它们，然后给我发电子邮件。

张量流概率

我们将使用 TensorFlow 2 作为基础库，它有一些类似于 PyTorch 的实现，尽管受益于我们将在后面章节中使用的许多特性。TensorFlow 通常也受益于语言和平台的可移植性。您可以从从 Java 到 C#的几乎任何语言中调用用 Python TensorFlow 编写的模型，并非常容易地将模型部署到移动设备上。

TensorFlow 和 TensorFlow-Probability 很容易安装。确保你已经安装了 Python(Python 3，如果你还没有这个或者任何 Python 发行版，我强烈推荐你从:www.anaconda.com下载 Anaconda 发行版，它包含了你在这里开始所需要的一切。)

安装张量流和张量流-概率:

如果要使用 Anaconda/Miniconda 安装 TensorFlow，请使用 conda 软件包管理器并键入:

conda install tensorflow

或者对于 GPU 支持:

conda install tensorflow-gpu

在撰写本文时，如果你想获得 GPU 支持，用 conda 而不是 pip 安装有一些好处。conda 版本自动安装数学库，显著提高 CPU 训练性能。使用 conda 安装 GPU(tensor flow-GPU)也更容易，因为默认情况下，它会安装 CUDA 和 cuDNN，而这两个软件需要使用 pip 进行独立安装。GPU versoins 要求机器上有 CUDA(即 NVIDIA) GPU。然而，conda 并不像 pip 那样包含那么多不同的库，当同时使用 pip 时很容易崩溃。

要使用 pip 安装计算(非 GPU)版本:

pip install — upgrade tensorflow==2.*

或者，如果您想使用带画中画的 GPU 版本:

pip install — upgrade tensorflow-gpu==2.*)

如果你不确定，选择第一个。2 号。*部分确保安装最新的 2 个版本。该系列大部分可以通过 TensorFlow 1.13 以上的版本完成，但需要在导入 tensor flow 后使用:tf.enable_eager_execution()激活急切执行(在 TF2 默认激活)，此外还需要在该系列后面的章节中进行一些其他更改。

注意: 在安装 TensorFlow-Probability 之前必须安装核心 TensorFlow 库，意外的是 TensorFlow-Probability 安装不会为你安装核心库作为依赖。

通过 conda 安装 TFP，包括:

conda install -c conda-forge tensorflow-probability

对于 pip:

pip install — upgrade tensorflow-probability

然后检查一切是否正常，创建一个新的 Jupyter 笔记本，IPython 实例或 Python 脚本，并添加:

import tensorflow as tf
import tensorflow_probability as tfp
dist = tfp.distributionsrv_normal = dist.Normal(loc=0., scale=3.)
print(rv_normal.sample([1]))

你会认出这一章开头的最后两行。你在对一个随机变量进行采样，所以作为输出，你应该得到一个张量，它包含了每次执行最后一行时不同的随机值。

安装 Pytorch 和 Pyro:

对于 Pyro，你首先需要安装 PyTorch。通常，这可以通过 conda 以下列方式最容易地完成:

 conda install pytorch-cpu torchvision-cpu -c pytorch

或者如果有 NIVIDA GPU 可供您使用:

 conda install pytorch torchvision cudatoolkit=9.0 -c pytorch

Pip 安装也是可能的，但不幸的是平台和的 CUDA 版本有很大不同，因此如果你只想使用 pip，建议你参考 pytorch.org 的安装指南。

Pyro 随后安装有:

conda install -c gwerbin pyro-ppl

或者

 pip install pyro-ppl 

**注意:**软件包‘Pyro’(没有-ppl)是完全不同的软件，与 Pytorch-Pyro 无关。

要开始使用 Pyro，创建一个 Jupyter 笔记本，Ipython 控制台或 python 脚本并键入:

import torch
import pyro
import pyro.distributions as dist
rv_normal = dist.Normal(loc=0., scale=3.)
print(rv_normal.sample([1]))

你会认出这一章开头的最后两行。你在对一个随机变量进行采样，所以作为输出，你应该得到一个张量，它包含了每次执行最后一行时不同的随机值。

摘要

你已经学习了随机变量和分布之间的关系，并看到了它的作用。您实现了一个位于 0、标度为 3 的正态分布。位置是分布平均值的另一个词，而标度是标准差的另一个词。然而，我们将坚持使用术语位置和比例，因为这些术语有助于我们在处理其他(非正态)分布时保持一致。

下周我们将学习如何实现贝叶斯神经网络，并进一步了解这些分布和随机变量的适用范围。我打赌你等不及了！订阅博客以获得关于下一版本的通知。

[1]:你可能也听说过 PyMC3，这是另一个很好的库，你可以用它来构建概率模型。然而，PyMC3 总体上特别倾向于统计建模，对使用神经网络的人不太友好(例如，自动签名和 GPU 支持更具挑战性)。

[2]:例如，TensorFlow 2 删除了设置图形所需的“会话”上下文，转而使用函数。

[3]:如果你不能选择使用哪个库，你会发现 TensorFlow-Probability 比 Pyro 更容易使用和理解。然而，这份关于 Pyro 的文档非常好，但是从神经网络的角度来解释 TFP 却很简单。

贝叶斯神经网络:TensorFlow 和 Pytorch 中的 2 个完全连接

原文：https://towardsdatascience.com/bayesian-neural-networks-2-fully-connected-in-tensorflow-and-pytorch-7bf65fb4697?source=collection_archive---------10-----------------------

贝叶斯神经网络

Photo by 青 晨 on Unsplash

本章继续贝叶斯深度学习系列。可用上一章 此处 和下一章 此处 **。**本章我们将探索传统密集神经网络的替代解决方案。这些备选方案将调用神经网络中每个权重的概率分布，从而产生单个模型，该模型有效地包含基于相同数据训练的神经网络的无限集合。我们将用这些知识来解决我们这个时代的一个重要问题:煮一个鸡蛋需要多长时间。

章节目标:

• 熟悉密集贝叶斯模型的变分推理
• 了解如何将正常的全连接(密集)神经网络转换为贝叶斯神经网络
• 了解当前实施的优点和缺点

这些数据来自煮鸡蛋的实验。煮鸡蛋的时间和鸡蛋的重量(以克为单位)以及切开鸡蛋后的发现一起提供。调查结果分为三类:欠熟、半熟和全熟。我们想从鸡蛋的重量和煮的时间得到鸡蛋的结果。这个问题非常简单，以至于数据几乎是线性可分的⁠.但也不尽然，因为没有提供鸡蛋的预煮寿命(冰箱温度或室温下的橱柜储存),你会看到这改变了烹饪时间。没有缺失的数据，我们无法确定打开鸡蛋时会发现什么。知道了我们有多确定我们可以影响这里的结果，就像我们可以影响大多数问题一样。在这种情况下，如果相对有把握鸡蛋没有煮熟，我们会在敲开之前多煮一会儿。

让我们先看看数据，看看我们在处理什么。如果你想亲自感受这种不同，你可以在github . com/doctor loop/Bayesian deep learning/blob/master/egg _ times . CSV获取数据。您将需要 Pandas 和 Matplotlib 来探索数据。(pip install —升级 pandas matplotlib)将数据集下载到您正在工作的目录中。如果不确定，可以从 Jupyter 笔记本的单元格中键入 pwd 来查找目录的位置。

https://gist . github . com/doctor loop/5a 8633691 f 912d 403 e 04 a 663 Fe 02 E6 aa

https://gist . github . com/doctor loop/21e 30 BD F16 D1 f 88830666793 f 0080 c 63

图 2.01 鸡蛋结果散点图

让我们看一下柱状图。

https://gist . github . com/doctor loop/2 a5e 95 a68a 29315 f 167 E0 e 875 E7 FAE 16

图 2.02 按结果划分的排卵次数直方图

我似乎不太擅长做我喜欢的半熟鸡蛋，所以我们看到了一个相当大的阶级不平衡，相对于半熟的可爱鸡蛋，半熟鸡蛋的数量是半熟鸡蛋的两倍，熟鸡蛋的数量是半熟鸡蛋的三倍。这种类别不平衡可能会给传统的神经网络带来麻烦，导致它们表现不佳，并且不平衡的类别大小是一个常见的发现。

注意，我们没有将 density 设置为 True (False 是默认值，因此不需要指定)，因为我们对比较实际数字感兴趣。而如果我们要比较从三个随机变量中的一个中采样的概率，我们会希望设置 density=True，以将数据总和的直方图归一化为 1.0。

直方图对于数据的离散化可能很麻烦，因为我们需要显式指定(或使用算法为我们指定)收集值的箱数。箱的大小极大地影响了数据的显示方式。作为一种替代方法，我们可以选择使用核密度估计，但在比较组时通常会更好，因为我们在这里使用的是小提琴图。violin 图是一种混合图，带有显示平滑分布的箱形图，便于比较。在 python 中，用 seaborn 库来绘制 violin 图既容易又漂亮(pip install —升级 Seaborn)。但是为了保持一致，我将在 matplotlib 中显示它。

https://gist . github . com/doctor loop/C5 BAC 12d 7 b 7 EBD 7758 cc 83 f 0 Fe 931 fc 0

图 2.03 按结局划分的鸡蛋次数小提琴图。中间的条形表示平均值。

好极了。现在让我们强调一下我们实际上要学习的神经网络之间的架构差异。我们将首先实现一个经典的密集架构，然后将其转换成一个强大的 Bayesianesque 式的对等物。我们使用基本库(TensorFlow 或 Pytorch)实现密集模型，然后使用附加库(TensorFlow-Probability 或 Pyro)创建贝叶斯版本。不幸的是，TensorFlow 实现密集神经网络的代码与 Pytorch 的代码非常不同，所以请转到您想要使用的库部分。

张量流**/张量流-概率**

源代码可从以下网址获得

【github.com/DoctorLoop/BayesianDeepLearning/ 号

https://gist . github . com/doctor loop/3a 67 FD 1197 a 355 f 68 f 45076 be 0074844

我们使用 Keras 来实现，因为这是我们后面的 TFP 模型所需要的。我们通过调用编译后的模型开始训练，如下所示。

https://gist . github . com/doctor loop/f6c 74509046068 AE 7 AC 37 EFE 00d 08545

你会看到 TensorBoard 是用来记录我们的训练的——利用 TensorBoard 作为其强大的模型调试手段和 TensorFlow 的关键功能之一是很好的。如果我们想在训练期间在单元格下输出训练进度，我们可以使用 Keras 输出，将 model.fit 的 verbose 参数设置为 1。Keras 输出的问题是它会打印每个时期的损失。这里我们训练了 800 个纪元，所以我们会得到一个又长又乱的输出。因此，最好创建一个定制的日志记录器来控制打印损失的频率。自定义记录器的源代码非常简单，可以在 Github 的笔记本中找到。TensorBoard 回调(和任何自定义记录器)都被传递给 callbacks 参数，但如果您不想只记录笔记本中的每个纪元，可以忽略这一点。Keras 的伟大之处在于，它会为我们拆分数据进行验证，从而节省时间。这里 10%的数据用于验证。有些人使用 25%的数据进行验证，但由于数据是神经网络最重要的考虑因素，数据集往往相当大，因此 10%的数据很好，并使模型更有可能达到我们的训练目标。Github 笔记本还展示了如何使用类权重来解决我们之前讨论的类不平衡问题。最后，如果培训时间很重要，请确保验证不频繁，因为验证比培训本身慢得多(特别是使用自定义记录器时，大约慢 10 倍)，从而使培训时间更长。

最终损耗在 0.15 左右，精度 0.85。如果你在 1200 个历元之后训练超过 800 个历元，你会发现验证的准确性在下降，尽管损失仍在下降。这是因为我们已经开始过度适应训练数据。当我们过度拟合时，我们使我们的模型依赖于特定的训练实例，这意味着它不能很好地概括我们看不见的验证数据或模型⁠.的预期真实世界应用过度拟合是传统神经网络的祸根，传统神经网络通常需要大型数据集和早期停止方法来缓解。

但我们是来学习解决方案而不是问题的！进入我们的贝叶斯神经网络阶段，它可以抵抗过度拟合，不会特别担心类别不平衡，最重要的是 在很少数据的情况下表现得非常好 。

让我们看一下模型的贝叶斯实现。

https://gist . github . com/doctor loop/f 00 c 18 f 591553685 e 06 a 79 fdb 4 b 68 e 0

它与我们的传统模型非常相似，只是我们使用了 TensorFlow-Probability 的 flipout 层。我们指定一个核散度函数，它是前面提到的 Kullback-Leibler 散度。大概就是这样。

Pytorch/Pyro

当比较传统的密集模型和贝叶斯等价模型时，Pyro 做的事情不同。对于 Pyro，我们总是先创建一个传统的模型，然后通过添加两个新的函数来进行转换。传统模型需要提供一种从权重分布中自动采样值的方法。采样值被插入到传统模型上的相应位置，以产生训练进度的估计。从上一章你会记得训练如何调节(训练)体重分布。在清晰的语言中，这意味着训练改变了正态分布，并改变了它们的比例来代表每个权重。然而，为了进行预测，我们仍然插入单一重量的固体值。虽然我们不能使用整个分布来进行预测，但是我们可以进行许多预测，每个预测使用不同的采样权重值来近似分布。这就是我们首先构建的密集模型非常适合的地方。

我们将从一个类开始来封装我们的密集模型。完整的源代码可以在 https://github.com/DoctorLoop/BayesianDeepLearning 的在线笔记本上找到。

https://gist . github . com/doctor loop/33b 192030388 c 46 d0a 2591459 b2f 6623

我们设置了三层，只有前两层使用激活功能。然后，我们指定了损失函数和优化器。请注意，对于这个模型，我们使用的是 torch.optim.Adam (Pytorch 优化器),而不是 Pyro 优化器。如果你试图在这里使用 Pyro，它会抛出错误，因为参数是不同的。如果你以前用过 Pytorch，我们的训练循环没有什么特别的。800 个纪元完全是多余的，但那又怎样。

https://gist . github . com/doctor loop/75854 f 26 bee 106 d9 b 596 B2 ee 0544 C5 D5

[Out]:
…
Test accuracy: 0.88
Final loss:0.173 at epoch: 799 and learning rate: 0.001

这里的表现非常好，速度也很快。如果你是 Pytorch 的新手，并且在其他地方有经验，你可能会对 softmax 的缺乏感到惊讶。在 Pytorch 中，需要注意的是 softmax 内置于 CrossEntropyLoss 中。

让我们继续升级它。这是我提到的两个新功能，模型和向导。

https://gist . github . com/doctor loop/a 4d cf 6 b 0 f 7 f1 ECE 43087534 ebdfcaa 08

我们将在本系列的另一篇文章中详细讨论这些函数的作用。简而言之，该模型明确声明了用于每个图层的分布，以替换点值。而指南声明了用于调节(训练)这些分布的变量。您会注意到这些函数看起来很相似，但仔细观察后，模型列出了权重和偏差的单独分布，然后指南列出了模型 中每个权重和偏差分布的平均值和标准差 的分布。这听起来有点玄乎——因为确实如此。把思考留到以后，现在就感受一下训练。

https://gist . github . com/doctor loop/5a 3837 c 384181 b 391 df 5927 BD D5 e 2 ab 5

训练这个美丽的怪物很快。我们得到的精度略高于密集模型 0.88–0.90，损失为 0.18。我在这里实现了一些额外的函数来使训练更容易(trainloader 和 predict 函数),并在 GitHub 笔记本中提供了完整的源代码。

贝叶斯深度学习的损失函数

我们不能通过随机变量反向传播(因为根据定义它们是随机的)。因此，我们欺骗和重新参数化这些分布，并让训练更新分布参数。随着这一根本性的变化，我们需要一种不同的方法来计算培训损失。

你可能已经知道传统模型中使用的负对数似然。它反映了模型生成数据的概率。现在不要担心这意味着什么。相反，我们通过取大量样本的平均值来估算这一损失。但是，随着我们的负对数可能性，我们结合了一个新的损失，利用计算中的分布。新的损失是 kull back-lei bler 散度(KL 散度),并提供了两个分布如何彼此不同的度量。我们会经常提到 KL 散度，并发现它在其他领域也很有用，包括 uncertainty⁴.指标负对数似然被加到 KL 散度上以得到 ELBO 损失(边际似然的预期下限，也称为变分自由能)。ELBO 损失允许我们在训练期间近似分布，并且在可训练性和训练时间方面受益巨大。

最后，作为总结，让我们来看看当我们做预测时会发生什么，毕竟这是我们在这里的目的。

https://gist . github . com/doctor loop/41 c 31 BF 30339 DC 1 B4 B4 c 7 e 69 a 001 b 9 BF

我们制作 5000 个样本(这可能需要几分钟，取决于您的计算机)。使用传统的神经网络，所有这些预测都将是完全相同的，因此这将是一个毫无意义的努力。但是我们不再使用传统的神经网络。

https://gist . github . com/doctor loop/1 ae0d 213d 3875 DD 569d 582 c 810769 fc 7

**图 2.04:**3 个测试输入的 5000 个样本的直方图

预测函数从我们用 TFP 或 Pyro 训练的主模型中采样多个不同的模型版本。我们输入的三个测试实例都产生了定义明确的峰值。第一个测试实例(红色直方图)的重量为 61.2 克，沸腾时间为 4.8 分钟。大多数时候，我们可以看到我们的模型预测它将是一个半熟的鸡蛋，但有 5%的时间它预测一个半熟的鸡蛋，有 2%的时间它认为它将是一个半熟的鸡蛋。该模型对此预测相当有信心，但不如第二个测试示例(绿色直方图)有信心，第二个测试示例几乎总是预测煮了 6 分钟的 53g 鸡蛋的煮熟结果。如果每次预测都不一样，我们如何从模型中获得一致性？我们只是取平均值。如果预测太易变，即每类中有三分之一，我们就不想做预测，除非我们可以根据信息采取行动，告诉用户模型对结果不确定(或者问问我们自己是否需要对模型/数据做更多的工作！)在这两种情况下，我们都获得了传统模型无法提供的非常有力的信息。因此，我们可以认为传统模型总是傲慢的，而我们的新模型是恰当的。很谨慎，但是自信正确的时候很自信。我们不需要成为计算机科学心理学家来欣赏更好的模型——人格。

总结

在这一章中，我们已经探索了基本贝叶斯模型的变化，并看到了伴随它的一些主要优势。我们已经在 TensorFlow-Probability 和 Pyro 中看到了如何做到这一点。虽然该模型功能齐全，但在现阶段它并不完美，也不是真正的贝叶斯模型。“更真实”的模型是否重要取决于你的环境。在随后的文章中，我们将讨论 softmax 的使用等缺陷，以及如何解决这些问题。在下一篇文章中，我们的主要关注点将是使用贝叶斯卷积深度学习进行图像预测。如果我的写作看起来更正式一点，我可能会更有信心在那里看到你！

1 因此，当更简单的模型可以做得很好的时候，用神经网络来解决问题有点滥用。

2 不确定这里的实际应用是什么，早餐咖啡馆？

3 如果你的想法是:哇，张量流概率代码看起来简单多了，因为它现在为你做了运算。也就是说，我们可以通过使用 Pyro autoguide 来简化 Pyro 代码。正如它所暗示的那样，这对你起着向导的作用。但我们在这里是为了学习贝叶斯深度学习的神奇之处，所以我们需要在某些时候曝光！

4 有趣的是，我们可以通过直接使用 KL 散度来制作具有合理性能的简单分类器，即，可以使用局部二进制模式算法和 KL 来构建图像模式分类器。在这种伙伴关系中，KL 比较输入图像和具有已知图案的图像之间的角和边缘的分布。

贝叶斯神经网络:3 贝叶斯 CNN

原文：https://towardsdatascience.com/bayesian-neural-networks-3-bayesian-cnn-6ecd842eeff3?source=collection_archive---------16-----------------------

贝叶斯神经网络

知道什么时候被骗的超级深度学习

图像由 DoctorLoop 提供

这是贝叶斯深度学习系列的第三章。上一篇可用 此处可用 。

我们已经知道神经网络是傲慢的。但是标准神经网络的另一个缺点是容易被欺骗。想象一下，一个 CNN 肩负着像人脸识别这样道德上有问题的任务。它从安全摄像头获取信息。但是，当安全摄像头发送一只狗、一个人体模型或一个儿童玩偶的照片时，会发生什么呢？CNN 仍然会输出被类似人脸的东西欺骗的分类。CNN 呼吁贝叶斯治疗，因为我们不希望我们的工作被愚蠢的错误所破坏，因为在错误分类的后果很严重的地方，我们想知道网络有多可靠。一个娃娃，或者任何一个真人，成为通缉犯的可能性有多大？我们需要知道。

章节目标

• 了解如何实现贝叶斯卷积模型
• 了解我们如何在从未见过的情况下识别坏的输入数据
• 理解贝叶斯神经网络的参数问题如何影响训练

正如我们在以前的文章中发现的，贝叶斯分析处理的是分布而不是单个值。我们在正态分布中看到了这一点，在正态分布中，我们得到了一个连续的浮点值，并返回了最有可能的平均值。分类的分布变成了离散的(钢琴键而不是小提琴弦)。对于概率，我们将得到一个特定的结果，比如一个类、一个指数(或者一个音符，如果我们挤压音乐类比的话)。这种分布是由我们模型中的逻辑信息决定的。

能够用少得多的数据表现得很好，同时具有更好的概括能力，这使得贝叶斯神经网络更受欢迎。而且这还没有考虑其他优势。然而，贝叶斯实现的一个缺点在本章变得很重要。贝叶斯实现需要更多的参数。考虑到整个分布替换了每个权重值参数，令人惊讶的是只需要两倍的参数。很容易看出相当具体的“两倍”数字来自哪里，因为重量分布大多是正态分布，而这些正态分布各自有两个参数。

在《离经叛道庄园》中，作者自然地创造了一个奇怪的新数据集来突出这个问题。我们发现自己在帮助父母解决一个重要的问题。父母总是对测量他们孩子的身高感兴趣。嗯，他们再也不用担心整天带着尺子到处跑了。我们将创建一个从图片中估计高度的模型。一个由 836 个婴儿和初学走路的孩子的剪影以及他们的身高组成的数据集已经被创建。因此，我们解决的不是分类问题，而是回归问题。我们的目标是返回一个与照片中孩子的身高相对应的浮点值。这是一个比上一章的分类练习稍微难一点的问题，并且由于蜘蛛的偶尔出现而变得更加困难。虽然训练集只包括有效的人类轮廓，但训练后我们会加入蜘蛛只是为了好玩。这就是有趣的地方，因为我们想避免返回昆虫的身高测量值，但是我们不允许有任何训练中的昆虫的照片。

用于训练的儿童原始剪影示例(左右两幅图像),位于不可用于训练的多毛蜘蛛剪影的两侧。这些剪影都被随机缩放至原始尺寸的 1/2–2 倍。

当然，蜘蛛通常比孩子小。因此，为了防止简单的基于大小差异的歧视，蜘蛛被升级以占据与儿童相当的空间。此外，为了让这项任务变得相当困难，孩子们被随机地重新调整了比例，所以他们的轮廓在原来的 1/2 到 2 倍之间。当然，当我们对孩子的身高感兴趣时，重新调整没有多大意义！但是这个精心设计的愚蠢问题完美地展示了贝叶斯深度学习的力量。您将看到这些模型如何很好地概括真实世界的情况，只有很少的训练示例，并且没有任何他们可能收到的损坏数据的示例！

让我们陷入张量流概率的问题。完整的代码和数据可以在 Jupyter 笔记本上在线获得:https://github.com/DoctorLoop/BayesianDeepLearning。首先，我们将定义架构。

贝叶斯卷积架构(https://gist . github . com/doctor loop/293 AE 5c C3 BDA a2 CCC 333d 9 b 216 eacc 301)

在第一行中，我们清除任何可能已经在内存中的会话，清空任何参数存储和变量，这样就不会干扰新的运行。接下来，我们定义一个 lambda 函数，帮助我们通过上一章讨论的 Kullback-Leibler (KL)散度更新损失。然后，我们将这个λ传递给每个卷积层，以便可以参考近似分布和我们的先验分布之间的差异来更新损耗。严格来说，这并不是绝对必要的，因为该层的默认参数几乎相同。然而，不同之处在于，默认参数仅获得 KL 散度，我们更进一步，将其除以示例总数(836)。默认实现将 epoch 的总 KL 应用于每个示例。但是我们更喜欢的是将总的历元 KL 的一部分应用到每个例子中，而不是每次都应用总的 KL。虽然两者都将训练，我们看到通过缩放损失更好的结果。自己实验看看。

实际模型的定义与任何其他 keras 序列一样。当然，我们使用的是卷积 2DFlipout 层(我们将在后面讨论)，而不是通常的 Conv2D。你可能会感到惊讶，我们只使用两个卷积层，而使用数百个卷积层已经足够了。我们使用两个只是因为结果令人印象深刻，对于这个问题，我们真的不需要更多。我们还在神经元层之间加入了两个最大池层，两者都有相当大的步幅/池大小。如果你有一个问题，需要特别敏感的像素完美的测量，你可能要尝试删除这些。当然，这样做的代价是硬件需求的增加，所以建议对两者进行比较。

最后一层是单个密集(贝叶斯)神经元，因为我们只对一个输出感兴趣。这个输出将是我们的测量。就这么简单。

最后用均方误差损失(MSE)编译模型。这是欺骗性的，因为虽然我们只指定了 MSE，但我们也在每一层上添加了 KL。然而，我们自己定义了吉隆坡，因为我们是独立的 Bayesianists 谁想要给 Keras 一个当之无愧的休息。当我们在培训期间打印损失时，我们将看到 KL 参与的证据。它明显不同于(大于)单独的 MSE。两者的区别就是我们的 KL。

培养

让我们开始训练，看看这次损失:

贝叶斯卷积模型的训练说明(https://gist . github . com/doctor loop/4b 10 c 410 a 709 e 0 DFD 71 ace 8 b 004255 BC

[Out]:
....
Epoch 250/250
151/151 [==============================] - 1s 4ms/step - loss: 12.5878 - mse: 5.1539 - val_loss: 16.3906 - val_mse: 8.9721

这里没有什么需要注意的。损耗比较高而批量比较低！

为了首先解决损失问题，我们会反复发现，与传统模型相比，贝叶斯模型的损失值是一个更差的模型性能指标。第一个原因是因为我们合并了至少两个损失。当然，我们感兴趣的是损失的变化，而不是明确的值，但即使这样，变化也不总是清晰的，因为我们经常在训练期间逐步改变两个损失的相对影响。我们将在后面的章节中讨论这些考虑因素。现在只需要记住，看到一个分类模型损失几千(！)同时拥有完善的验证指标。

虽然有些人可能会嘲笑我微不足道的批量，并认为资源是稀缺的——他们错得不能再错了。使用贝叶斯模型，批量大小对训练的影响比我们预期的要大得多。这是我们经常认为我们理解的神经网络理论的一些领域的一个例子，但这将要求我们重新审视我们的信念。我们通常认为批量大小对训练速度至关重要。一些人也欣赏更大批量带来的减少的差异。然而使用贝叶斯模型，批量大小直接影响培训绩效。通过重复运行批量为 5 和 50 的相同模型来观察一下。您会注意到，当批量大小为 50 时，epochs 当然会快得多，但我们从未获得像批量大小为 5 时那样好的损失或性能指标。这不是一个小差异，而是巨大的差异！这很重要，因为很快就会清楚，批量大小是贝叶斯深度学习成功的一个非常有影响力的超参数。

很快就会清楚，批量大小是贝叶斯深度学习成功的一个非常有影响力的超参数

虽然起初我们有另一个超参数要优化，这似乎令人沮丧，但我们会发现自己能够通过一个非常简单的架构变化来提升性能，这比我们过去依赖的架构更简单(在本文底部的附录中，我们将讨论像 Flipout 这样推动这些变化的层)。

推理

最后我们得出结论。我们对从贝叶斯主模型中做出多重预测感兴趣。每个输出都略有不同，因为每个预测都将由一个新的模型做出，该模型填充了从我们训练的贝叶斯主模型的权重分布中采样的权重。

两个列表理解，每个为两个不同的输入图像生成 1000 个预测。https://gist . github . com/doctor loop/09552736976 a7e 0 a 32 E3 f 27d 28 a4 ee1 c

在上面的代码中，我们使用列表理解方式进行循环预测。如果我们只提供一个输入数组(1000 x 126 x 126 x 1)并一次做出所有的预测，不是更快吗？事实上这会快得多。但与此同时，它会违背目的，因为它是一个单独的模型。predict 调用从我们的贝叶斯训练的分布中采样新的权重。因此，每个预测调用负责创建一个独特的新模型，该模型受我们在培训中创建的分布的约束。如果我们输入 1000 张图片，只进行一次预测调用，所有的预测都将是相同的，因为我们将使用单个权重样本，从而模拟一个标准模型。我们更感兴趣的是利用贝叶斯训练创建的无限模型包的能力。我们称之为模型集合，我们利用多个不同模型的集合来获得相同数据的许多不同观点。许多观点的一致是最重要的，它告诉我们输入数据的质量。

根据有效输入(baby_predictions)和无效输入(spider_predictions)绘制预测。https://gist . github . com/doctor loop/41 Abe 385934 fb0f 7728 ba 048564 e26d 4

针对单个有效输入(绿色)的 1000 个高度预测和单个无效输入(红色)的 1000 个高度预测的密度图。当模型被给予具有蜘蛛的无效输入时，预测的分布显示了预测之间的不一致，表明高度不确定性。对有效婴儿轮廓的测量预测的相似性显示了指示有把握的预测的预测一致性。

在上面的代码和图中，我们制作了一张婴儿图片(绿色)和一张蜘蛛图片(红色)的 1000 个高度预测的密度图。我们可以看到，对婴儿身高的预测非常紧密地聚集在 51 像素左右(其平均值和期望值)。而大约 30%的预测正好是这个测量值(恰好是真实值),并且几乎所有的预测都在事实的一个像素之内！另一方面，虽然对蜘蛛的预测也集中在一个值(90 像素)上，但只有不到 4%的预测是在期望值上，并且预测在从 51 像素一直到 134 像素的范围内分散得更广。显然蜘蛛上的预言彼此不一致。因此，我们可以凭直觉认为，我们的贝叶斯模型对于无效对象的预测是不确定的，而我们的贝叶斯模型对于来自训练的对象的相关预测是有信心的。这正是我们想要的。

在下一篇文章中，我们将探索如何使简单的贝叶斯模型比复杂的标准模型更好。我们还将了解如何利用其他类型的不确定性来指导培训，以及如何优化和比较模型以找到最佳方案。

附录:张量流-概率卷积层

如果你最近阅读了文档或任何论文，你可能已经找到了解决贝叶斯深度学习的不同方法。TensorFlow Probability 为卷积层实现了两种方法(更多方法可用于密集层)。第一个是重新参数化层:

TFP . layers . convolution 2d 重新参数化。重新参数化允许我们通过分布的最可能值来计算梯度。因此，我们操纵描述分布的参数，而不是神经网络中的权重值。处理分布参数意味着实际的分布可以被忽略，并被有效地抽象掉。描述分布的参数可以被认为是分布对象的替身，就像纸币代表黄金等实物资产一样。在这两种情况下，替身都是首选，因为这样更方便。在训练中，我们方便地避免了通过随机变量尝试反向传播的尴尬(尴尬是因为它不起作用)。

重新参数化速度很快，但遗憾的是，实际上需要将一批实例的所有权重设置为相同的值。如果重量被单独记录而不是共享，内存需求将会激增。共享权重是有效的，但是增加了方差，使得训练需要更多的历元。

flipout 层:TFP . layers . convolution 2d flipout 采用不同的方法。虽然类似，但它受益于一个特殊的损耗梯度估计器。这种 flipout 估计器在小批量中摇动权重，使它们更加相互独立。反过来，与重新参数化方法相比，抖动减少了方差，并且需要更少的训练时期。但是有一个问题。虽然 flipout 需要更少的历元，但它实际上需要两倍的计算！幸运的是，这些计算可以并行化，但我们仍然会发现一个模型每个时期需要 25–50%的时间(取决于硬件)，即使训练总共需要更少的时期。

在没有重新参数化分布的情况下，我们打破了大样本给出更好估计的假设。虽然我们中的许多人并没有从这些方面考虑训练，但我们一直依赖于这个假设。因此，通过重新参数化，我们描述了最可能值的变化，而不是样本中最可能的变化，我们无法预测样本中最可能的变化，因为变量不是随机的，如果我们可以预测的话。

具有张量流概率的贝叶斯神经网络

原文：https://towardsdatascience.com/bayesian-neural-networks-with-tensorflow-probability-fbce27d6ef6?source=collection_archive---------2-----------------------

概率模型不确定性预测的逐步指南。

图片 via 维基百科。

介绍

机器学习 模型通常从数据开发为 确定性 机器，其使用通过 最大似然 方法计算的参数权重的点估计来映射输入到输出。然而，在后台有很多统计上的侥幸。例如，数据集本身是来自叠加了加性噪声的未知分布的任意大小的点的有限随机集合，并且对于这种特定的点集合，不同的模型(即不同的参数组合)可能是合理的。因此，对于正在进行的参数和预测，存在一些 不确定性 。 贝叶斯统计 提供了一个框架来处理所谓的不确定性，随着 张量流概率 、 概率建模 的发布，这已经变得容易多了，我将在这篇文章中演示。要知道不会提供任何理论背景；关于这个话题的理论，我真的可以推荐 Gelman 等人的书《贝叶斯数据分析》,这本书是免费的 PDF 文件。

**** [## “贝叶斯数据分析”这本书的主页

这是本书的 pdf 格式，可以下载用于非商业目的。

www.stat.columbia.edu](http://www.stat.columbia.edu/~gelman/book/)

贝叶斯神经网络

贝叶斯神经网络 的特征在于其在权重(参数)和/或输出上的分布。取决于是否考虑了听觉的、认知的或两者都考虑的不确定性，贝叶斯神经网络的代码看起来略有不同。为了演示工作原理，来自 De Vito 的空气质量数据集将作为一个例子。它包含来自不同污染物化学传感器的数据(如电压)以及作为一年时间序列的参考数据，这些数据是在一个汽车交通繁忙的意大利城市的主要街道上收集的，目标是构建从传感器响应到参考浓度的映射(图 1)，即构建一个校准函数作为 回归 任务。

图 1: 参考和传感器数据。注意参考数据中的相关性。(图片由作者提供)

启动

如果您还没有安装 TensorFlow Probability，您可以使用 pip 来完成，但是在此之前创建一个虚拟环境可能是一个好主意。(由于命令在以后的版本中可能会改变，您可能希望安装我用过的那些。)

**# Install libraries.**pip install tensorflow==2.1.0
pip install tensorflow-probability==0.9.0

打开你最喜欢的编辑器或者 JupyterLab。导入所有必要的库。

**# Load libriaries and functions.**import pandas as pd
import numpy as np
import tensorflow as tf
tfk = tf.keras
tf.keras.backend.set_floatx("float64")
import tensorflow_probability as tfp
tfd = tfp.distributions
from sklearn.preprocessing import StandardScaler
from sklearn.ensemble import IsolationForest**# Define helper functions.**scaler = StandardScaler()
detector = IsolationForest(n_estimators=1000, behaviour="deprecated", contamination="auto", random_state=0)
neg_log_likelihood = lambda x, rv_x: -rv_x.log_prob(x)

接下来，获取数据集(链接可以在上面找到)并将其作为熊猫数据框架加载。由于传感器容易因老化而漂移，最好丢弃过去六个月的数据。

**# Load data and keep only first six months due to drift.**data = pd.read_excel("data.xlsx")
data = data[data["Date"] <= "2004-09-10"]

预处理

数据相当杂乱，必须先进行预处理。我们将关注大部分时间测量的输入和输出(一个传感器很早就死了)。删除缺少值的行后，对数据进行缩放。之后，使用 隔离林 检测并移除异常值。

**# Select columns and remove rows with missing values.**columns = ["PT08.S1(CO)", "PT08.S3(NOx)", "PT08.S4(NO2)", "PT08.S5(O3)", "T", "AH", "CO(GT)", "C6H6(GT)", "NOx(GT)", "NO2(GT)"]
data = data[columns].dropna(axis=0)**# Scale data to zero mean and unit variance.**X_t = scaler.fit_transform(data)**# Remove outliers.**is_inlier = detector.fit_predict(X_t)
X_t = X_t[(is_inlier > 0),:]**# Restore frame.**dataset = pd.DataFrame(X_t, columns=columns)**# Select labels for inputs and outputs.**inputs = ["PT08.S1(CO)", "PT08.S3(NOx)", "PT08.S4(NO2)", "PT08.S5(O3)", "T", "AH"]
outputs = ["CO(GT)", "C6H6(GT)", "NOx(GT)", "NO2(GT)"]

数据处理

TensorFlow 提供了一个 数据集类 来构造训练集和测试集。我们将使用 70%的数据作为训练集。集合被混洗并且重复的批次被构造。

**# Define some hyperparameters.**n_epochs = 50
n_samples = dataset.shape[0]
n_batches = 10
batch_size = np.floor(n_samples/n_batches)
buffer_size = n_samples**# Define training and test data sizes.**n_train = int(0.7*dataset.shape[0])**# Define dataset instance.**data = tf.data.Dataset.from_tensor_slices((dataset[inputs].values, dataset[outputs].values))
data = data.shuffle(n_samples, reshuffle_each_iteration=True)**# Define train and test data instances.**data_train = data.take(n_train).batch(batch_size).repeat(n_epochs)
data_test = data.skip(n_train).batch(1)

模型结构

听觉不确定性

为了解决由输出中的噪声引起的听觉不确定性，密集图层与概率图层相结合。更具体地，输出的均值和协方差矩阵被建模为输入和参数权重的函数。第一个隐藏层应包括 10 个节点，第二个隐藏层需要 4 个节点用于均值，加上 10 个节点用于最终层中的四维(有 4 个输出)多元高斯后验概率分布的方差和协方差。这是通过使用最后一层 ( MultivariateNormalTriL)的params_size 方法来实现的，这是 后验概率分布 结构的声明，在这种情况下是多元正态分布，其中仅估计了协方差矩阵的一半(由于对称性)。模型中的参数总数为 224 个——由 变分法 估算。这个神经网络的确定性版本由一个输入层、十个 潜变量 (隐藏节点)和一个输出层(114 个参数)组成，其中不包括参数权重的不确定性。

**# Define prior for regularization.**prior = tfd.Independent(tfd.Normal(loc=tf.zeros(len(outputs), dtype=tf.float64), scale=1.0), reinterpreted_batch_ndims=1)**# Define model instance.**model = tfk.Sequential([
tfk.layers.InputLayer(input_shape=(len(inputs),), name="input"),
tfk.layers.Dense(10, activation="relu", name="dense_1"),
tfk.layers.Dense(tfp.layers.MultivariateNormalTriL.params_size(
len(outputs)), activation=None, name="distribution_weights"),
tfp.layers.MultivariateNormalTriL(len(outputs), activity_regularizer=tfp.layers.KLDivergenceRegularizer(prior, weight=1/n_batches), name="output")
], name="model")**# Compile model.**model.compile(optimizer="adam", loss=neg_log_likelihood)**# Run training session.**model.fit(data_train, epochs=n_epochs, validation_data=data_test, verbose=False)**# Describe model.**model.summary()

activity_regularizer 参数作为输出层的先验(权重必须根据批次数进行调整)。根据您机器的规格，培训可能需要一段时间。该算法需要大约 50 个历元才能收敛(图 2)。

图 2: 有训练和验证损失的学习曲线。(图片由作者提供)

听觉和认知的不确定性

考虑到听觉和认知的不确定性(参数权重的不确定性)，密集层必须与 翻转 层(DenseFlipout)或 变化 层(DenseVariational)交换。这种模型有更多的参数，因为每个权重都由具有非共享平均值和标准偏差的正态分布来参数化，因此参数权重的数量加倍。对于不同的预测，权重将被重新采样，在这种情况下，贝叶斯神经网络将表现得像一个 集合 。

tfp.layers.DenseFlipout(10, activation="relu", name="dense_1")

权重的默认先验分布是tfd.Normal(loc=0., scale=1.)，可以使用kernel_prior_fn 参数进行调整。

预言；预测；预告

因为是概率模型，所以执行 蒙特卡罗实验 来提供预测。特别是，样本 x 的每一次预测都会产生不同的输出 y ，这就是为什么必须计算许多单独预测的期望值。此外，可以通过这种方式确定方差。

**# Predict.**samples = 500
iterations = 10
test_iterator = tf.compat.v1.data.make_one_shot_iterator(data_test)
X_true, Y_true, Y_pred = np.empty(shape=(samples, len(inputs))), np.empty(shape=(samples, len(outputs))), np.empty(shape=(samples, len(outputs), iterations))
for i in range(samples):
    features, labels = test_iterator.get_next()
    X_true[i,:] = features
    Y_true[i,:] = labels.numpy()
    for k in range(iterations):
        Y_pred[i,:,k] = model.predict(features)

**# Calculate mean and standard deviation.**Y_pred_m = np.mean(Y_pred, axis=-1)
Y_pred_s = np.std(Y_pred, axis=-1)

图 3 显示了所有输出的测量数据与预测值的对比。 决定系数 约为 0.86，斜率为 0.84——不算太差。

图 3: 测量值与预测值的一致性。(图片由作者提供)

预测的不确定性可以通过绘制误差线和期望值来可视化(图 4)。在这种情况下，误差棒是标准偏差的 1.96 倍，即占概率的 95%。

图 4: 数据与预测的比较。(图片由作者提供)

理论上，由于额外的不确定性信息，贝叶斯方法优于确定性方法，但由于其高计算成本，并不总是可行的。最近的研究围绕着开发新的方法来克服这些限制。

[## 用贝叶斯推理量化预测的不确定性

贝叶斯统计最重要的特征之一。

towardsdatascience.com](/uncertainty-quantification-of-predictions-with-bayesian-inference-6192e31a9fa9)****