抛硬币，随机变量的算术，和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很快就变得有用了。也就是说，我们可能需要对随机变量本身进行操作，并且只有在计算图建立之后才开始采样。在lea中，正是这样做的（尽管只针对离散分布），上面的例子变成

^{pr2}$

看起来很有魅力。输入PyMC3，这是概率编程中比较流行的库之一。现在，PyMC3打算用于MCMC和Bayesian建模，但是它有我们上面实验所需的构建块。唉

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。不幸的是，我缺乏理论上的理解，为什么吉布斯采样器会产生这种特殊的结果（我真的应该去看看书）。使用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的模拟失败了，在普通的、非MC的MC中使用PyMC3有什么缺陷吗？我正在尝试做的甚至是可能的首先是PyMC3？在