抛硬币、随机变量的算术和 PyMC3答案

【问题标题】：Coin tosses, arithmetic of random variables, and PyMC3抛硬币、随机变量的算术和 PyMC3
【发布时间】：2018-07-27 19:57:45
【问题描述】：

我发现自己想在 Python 中执行随机变量的算术运算；为了举例，让我们考虑反复抛两个独立的公平硬币并计算正面数量的实验。

使用scipy.stats 可以直接从每个随机变量中独立采样，我们可以立即开始获取结果

In [5]: scipy.stats.bernoulli(0.5).rvs(10) + scipy.stats.bernoulli(0.5).rvs(10)
Out[5]: array([1, 0, 0, 0, 1, 1, 1, 2, 1, 2])

现在，悲观主义者会说我们甚至不必走那么远，而可以只做np.random.randint(2, size=10) + np.random.randint(2, size=10)，而愤世嫉俗的人会注意到我们可以只计算总和而不必采样任何东西。

他们是对的。所以，假设我们有更多的变量和更复杂的操作要对它们执行，graphical models 很快就会变得有用。也就是说，我们可能希望对随机变量本身进行操作，并且仅在设置计算图时才开始采样。在 lea 中，正是这样做的（尽管仅适用于离散分布），上面的示例变为

In [1]: from lea import Lea

In [7]: (Lea.bernoulli(0.5) + Lea.bernoulli(0.5)).random(10)
Out[7]: (0, 2, 0, 2, 0, 2, 1, 1, 1, 2)

看起来像一个魅力。输入PyMC3，它是概率编程最流行的库之一。现在，PyMC3 特别适用于 MCMC 和贝叶斯建模，但它具有我们上述实验所需的构建块。唉，

In [1]: import pymc3 as pm

In [2]: pm.__version__
Out[2]: '3.2'

In [3]: with pm.Model() as model:
   ...:     x = pm.Bernoulli('x', 0.5)
   ...:     y = pm.Bernoulli('y', 0.5)
   ...:     z = pm.Deterministic('z', x+y)
   ...:     trace = pm.sample(10)
   ...:
Assigned BinaryGibbsMetropolis to x
Assigned BinaryGibbsMetropolis to y
100%|███████████████████████████████████████| 510/510 [00:02<00:00, 254.22it/s]

In [4]: trace['z']
Out[4]: array([2, 0, 2, 0, 2, 0, 2, 0, 2, 0], dtype=int64)

Not exactly random。不幸的是，我缺乏对 Gibbs 采样器为什么会产生这种特殊结果的理论理解（实际上我可能应该直接阅读书籍）。使用 step=pm.Metropolis() 代替，我们在一天结束时得到正确的分布，即使个别样本与其邻居有很强的相关性（正如 MCMC 所预期的那样）。

In [8]: with pm.Model() as model:
   ...:     x = pm.Bernoulli('x', 0.5)
   ...:     y = pm.Bernoulli('y', 0.5)
   ...:     z = pm.Deterministic('z', x+y)
   ...:     trace = pm.sample(10000, step=pm.Metropolis())
   ...:
100%|██████████████████████████████████████████████████████████████████████████████████████████| 10500/10500 [00:02<00:00, 5161.18it/s]

In [14]: collections.Counter(trace['z'])
Out[14]: Counter({0: 2493, 1: 5024, 2: 2483})

所以，也许我可以继续使用pm.Metropolis 来模拟我的算术后分布，但我担心我遗漏了什么，所以问题最终变成了：为什么step-上面的模拟失败了，将 PyMC3 用于普通、非 MC、MC 是否有任何陷阱，我首先尝试在 PyMC3 中做些什么？

【问题讨论】：

这个我还没调试完，但初步看起来像BinaryGibbsMetropolis子类ArrayStep，似乎实现了多个随机变量一起更新的行为。我认为这意味着存在一个错误，并且在某处将变量x 和y 分配给BinaryGibbsMetropolis 的同一个实例，这就是它们匹配的原因在每个随机步骤中彼此（Gibbs 采样器更新整个变量向量，而不是单独更新组件或使用两个单独的采样器处理变量）。
即使我使用这样的东西，我仍然会看到同样奇怪的固定行为：pm.sample(10000, step=[pm.BinaryGibbsMetropolis([x]), pm.BinaryGibbsMetropolis([y])], random_seed=[109, 287])。但是，如果我什至只将其中一个更改为普通的BinaryMetropolis 或Metropolis，那么它可以解决问题。你怎么敢在星期五晚上这样对我！
@ely：感谢您对此进行调查。所以特别是你似乎在说第一次尝试应该真的奏效了？
是的，我认为它应该可以工作，而且我在BinaryGibbsMetropolis 和BinaryMetropolis 的astep 成员函数中看不到任何不同的代码，所以@987654345 @ 的接受率低于 100%，而 BinaryGibbsMetropolis 似乎有 100% 的接受率，从代码来看似乎是不可能的，所以肯定有其他事情发生。
不应该是pm.__version__吗？

标签： python random probability montecarlo pymc3

【解决方案1】：

colcarroll 的评论：

[2 月2018 年 2 月 21 日]：绝对是一个错误 - github.com/pymc-devs/pymc3/issues/2866。你正在做的应该工作，但不是图书馆的意图。您将使用 PyMC3 来推理不确定性（也许观察 z 并推理 x 和 y 的概率）。我认为您的前两种方法，也许石榴库会更有效。请参阅 stackoverflow.com/questions/46454814/... –

[2 月2018 年 2 月 25 日]：这已由 Junpeng Lao 修复在 master 上（参见 github.com/pymc-devs/pymc3/pull/2867）。请参阅 andrewgelman.com/2018/01/18/... 了解“反相关抽奖”的背景。我不确定 stackoverflow 想要如何处理这样的问题。

【讨论】：