Comparing two biased coins (newbie example from Krushke's book)

I am a complete newbie to Bayesian statistics and MCMC, so I worked my way through Bayesian Data Analysis: A Tutorial with R and BUGS by John Krushk. To test my understanding, I am trying to translate his examples from BUGS to PyMC.

In ch. 8, he observes a sequence of flips from each of two (potentially biased) coins and tries to estimate the difference in their offsets. Below is the bias for each coin theta

, and the observed flips y

:

import numpy as np
import pymc
priorA = 3
priorB = 3
theta1 = pymc.Beta('theta1', alpha=priorA, beta=priorB)
theta2 = pymc.Beta('theta2', alpha=priorA, beta=priorB)
y1 = pymc.Bernoulli('y1', p=theta1, value=np.array([1,1,1,1,1,0,0]), observed=True)
y2 = pymc.Bernoulli('y2', p=theta2, value=np.array([1,1,0,0,0,0,0]), observed=True)

      

        
        
        
      

    

Of course, the two coins are independent. If I simulate them separately, then look at the difference in thetas, I get the same answer as the book. (It is also consistent with analytic solution and grid integration.)

m1 = pymc.MCMC([y1, theta1])
m1.sample(iter=10000, burn=1000, thin=2)
m2 = pymc.MCMC([y2, theta2])
m2.sample(iter=10000, burn=1000, thin=2)
dt = m1.trace('theta1')[:] - m2.trace('theta2')[:]
print np.average(dt), pymc.utils.hpd(dt, alpha=0.05)
# 0.23098692189 [-0.15303663  0.55686801]

      

        
        
        
      

    

If, on the other hand, I try to simulate two coins at the same time, as part of the same model, I get a very different answer that does not agree well with any other method.

model = pymc.MCMC([y1, theta1, y2, theta2])
model.sample(iter=10000, burn=1000, thin=2)
dt = model.trace('theta1')[:] - model.trace('theta2')[:]
print np.average(dt), pymc.utils.hpd(dt, alpha=0.05)
# 0.330061191424 [-0.09999254  0.7190196 ]

      

        
        
        
      

    

So my first question is why? From what little I understand from theory, my wild guess is that MCMC is having difficulty finding a two-parameter model. (BUGS seems to handle this just fine, though.)

It's really weird that I was doing all this in iPython notebook and it looks like there was a bug in PyMC. If I run the independent coin model, reload my kernel (Kernel | Restart or File | Close and Halt) and then run the co-coin model, co-coins will give a response similar (but not identical) to independent coins (dtheta on average ~ 0.23). If I run the models in reverse order (restarting the kernel in between), they both produce a (wrong) response from the co-coin model, average dtheta ~ 0.33. I can only get two models to get different answers if I shut down my iPython notebook server completely in between. Since this will also unload all shared libraries from memory, I would expect this to mean that the Fortran / C part of PyMC caches some of these models in a memory location.which is shared with Python interpreter instances. Versions are Numpy 1.8.2, PyMC 2.3.3, Python 2.7.8, iPython 2.1.0, Anaconda 2.0.0.

I would really appreciate any help in understanding what's going on here. I understand that these are stupid, trivial models, but the strange behavior of PyMC does not inspire confidence at the moment!