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

【问题标题】：How to model sum of two Poisson distributions with pymc 2?如何使用 pymc 2 对两个泊松分布的总和进行建模？
【发布时间】：2015-02-13 05:26:39
【问题描述】：

我正在尝试使用 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 采样器，我得到：

使用粒子吉布斯，我得到：

 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。

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

【问题讨论】：

您能否更详细地描述英国国教模型的输出？我不熟悉这个系统。
是的，它基本上是频率计数。例如9977 (7 0) 表示 a=7, b=0 在 10000 个样本中出现了 9977 次。

标签： python probability pymc

【解决方案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) 的东西时，PyMC 做了一些过于聪明的事情
precision 是一个数字，不是随机数，是一个大数字，而不是一个小数字。这是因为 PyMC 使用精度来表示正态分布中的 1/方差，但在英国国教（我认为）中，精度意味着您需要等式运算符有多接近精确度。
mcmc 采用自适应 Metropolis step 方法，老化时间较长。这一点很重要，因为a 和b 的联合后验分布具有极强的相关性，除非它弄清楚这一点，否则 MCMC 步骤将不会去任何地方。

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

【讨论】：