如何使用 pymc 2 对两个泊松分布的总和进行建模？

Question

我正在尝试使用 pymc 2 对一个简单的概率编程示例进行建模。我一直在使用其他语言，例如 Church 和 Anglican，并且能够毫无困难地对这个问题进行建模。但是，我似乎无法在 Python.

中弄清楚

这是 code in Anglican，我认为它很容易理解：

[assume a (- (poisson 100) 100)]
[assume b (- (poisson 100) 100)]
[observe (normal (+ a b) .00001) 7]
[predict (list a b)]

使用 Metropolis-Hastings 采样器，我得到：

使用 Particle Gibbs，我得到：

 669 (-1 8)
  71 (-10 17)
  66 (-11 18)
 208 (-12 19)
  19 (-13 20)
  84 (-14 21)
  72 (-15 22)
 441 (-2 9)
...and so on...

我正尝试在 pymc 中对此进行建模，如下所示：

def make_model():
    a = (pymc.Poisson("a", 100) - 100)
    b = (pymc.Poisson("b", 100) - 100)

    precision = pymc.Uniform('precision', lower=.0001, upper=1.0)

    @pymc.deterministic
    def mu(a=a, b=b):
        return a+b

    y = pymc.Normal("y", mu=mu, tau=precision, observed=True, value=7)

    return pymc.Model(locals())

def run_mcmc(model):
    mcmc = pymc.MCMC(model)
    mcmc.sample(5000, burn=1000, thin=2)
    return mcmc

result = run_mcmc(make_model())
pymc.Matplot.plot(result)

我得到 a 和 b 大约为 100 的轨迹。但是，如果我运行 (pymc.Poisson("a", 100) - 100).value，我得到的数字更接近 0。

我是不是漏掉了什么？我对这些可能性感到很兴奋，但现在我很困惑！感谢您的帮助！

Answer 1

如果我理解正确的话，这是一个很好的例子来说明 Anglican 和 PyMC 在思维上的一些差异。

我认为这是您的 PyMC 代码的调整版本：

def make_model():
    a = pymc.Poisson("a", 100)  # better to have the stochastics themselves available
    b = pymc.Poisson("b", 100)

    precision = 1e-4**-2 #  Seems like precision is fixed in Anglican model (and different from the meaning of precision in PyMC)

    @pymc.deterministic
    def mu(a=a, b=b):
        return (a-100) + (b-100)

    y = pymc.Normal("y", mu=mu, tau=precision, observed=True, value=7)

    return pymc.Model(locals())

def run_mcmc(model):
    mcmc = pymc.MCMC(model)
    mcmc.use_step_method(pymc.AdaptiveMetropolis, [mcmc.a, mcmc.b])
    mcmc.sample(20000, burn=10000, thin=10)
    return mcmc

result = run_mcmc(make_model())
pymc.Matplot.plot(result)

以下是我的代码中的主要区别：

a 和 b 是随机指标。当你在你的版本

(stochastic - 100)

precision是一个数字，不是随机数，是一个大数字，不是一个小数字。这是因为 PyMC 使用精度表示正态分布中的 1/方差，但在 Anglican 中（我认为）精度表示您需要等式运算符的精度有多接近。
mcmc采用自适应Metropolis步法，老化时间较长。这很重要，因为 a 和 b 的联合后验分布具有极大的相关性，除非弄清楚这一点，否则 MCMC 步骤不会走到任何地方。

这里 an IPython Notebook 显示了更多细节。

如何使用 pymc 2 对两个泊松分布的总和进行建模？

How to model sum of two Poisson distributions with pymc 2?

python

probability

pymc