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

我发现自己想要在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 . 不幸的是,我对吉布斯取样器为什么会产生这种特殊结果缺乏理论上的理解(实际上我应该认真研究一下)。使用 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中首先做的是什么?