抛硬币,随机变量的算术,和PyMC3

2024-06-15 20:56:11 发布

您现在位置:Python中文网/ 问答频道 /正文

我发现自己想在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?在


Tags: inasstatssteptracerandom算术scipy