抛硬币、随机变量的算法和 PyMC3

Question

我发现自己想在 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 quickly become useful. That is, we might want to operate on the random variables themselves and only start sampling when our graph of computation is set up. In 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 来模拟我的 post-算术分布，但我担心我遗漏了什么，所以问题最终变成了：为什么上面的 step-less 模拟会失败，并且在将 PyMC3 用于普通的、非 MC、MC 时是否存在任何缺陷，而我首先尝试在 PyMC3 中做的事情甚至可能吗？

Answer 1

colcarroll的评论：

[Feb. 21, 2018]: Definitely a bug - github.com/pymc-devs/pymc3/issues/2866 . What you are doing should work, but is not the intention of the library. You would use PyMC3 to reason about uncertainty (perhaps observing z and reasoning about the probabilities of x and y). I think your first two approaches, and perhaps the pomegranate library would be more efficient. See whosebug.com/questions/46454814/… –

[Feb. 25, 2018]: This is now fixed on master (see github.com/pymc-devs/pymc3/pull/2867) by Junpeng Lao. See andrewgelman.com/2018/01/18/… for background on "Anticorrelated draws". I am not sure how Whosebug wants to handle a question like this.

抛硬币、随机变量的算法和 PyMC3

Coin tosses, arithmetic of random variables, and PyMC3

python

random

probability

montecarlo

pymc3