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

Question

I'm trying to model a simple probabilistic programming example using pymc 2. I've been playing with other languages such as Church and Anglican and am able to model this problem without difficulty. However, I can't seem to figure it out in Python.

Here is the code in Anglican, I think it's pretty self-explanatory:

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

Using the Metropolis-Hastings sampler, I get:

   1 (10 1)
   2 (10 8)
9977 (7 0)
  20 (7 1)

With Particle Gibbs, I get:

 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...

I'm trying to model this in pymc like so:

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)

I'm geting traces where a and b are around 100. However, if I run (pymc.Poisson("a", 100) - 100).value, I get numbers closer to 0.

Am I missing something here? I'm excited about the possibilities, but am very confused at the moment! Thanks for any help!

Answer 1

If I understand this correctly, it is a nice example to demonstrate some of the differences in thinking between Anglican and PyMC.

Here is a tweaked version of your PyMC code that I think captures your intention:

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)

Here are the key differences in my code:

• a and b are stochastics. PyMC does something too clever for its own good when you use the (stochastic - 100) stuff in your version
• precision is a number, not a stochastic, and a big number, not a little number. This is because PyMC uses precision to mean 1/variance in a normal distribution, but in Anglican (I think) precision means how close to precise do you need the equality operator to be.
• mcmc uses the adaptive Metropolis step method, with a longer burn-in period. This is important because the joint posterior distribution of a and b has extreme correlation, and the MCMC steps will not go anywhere unless it figures this out.

Here is an IPython Notebook which shows a bit more detail.