I recently completed a class that made heavy use of PyMC3 and have been thinking of ways to use MCMC sampling in my PhD research. I have been trying to use PyMC3 for fitting some experimental data that comes from single quantum emitters, which works well for the simplest model but I have been having issues when I increase the complexity slightly.
I have data representing samples taken at particular time points (t_vec), and I have a full set of these time points for each experiment. So I have a matrix of shape (N_EXP, len(time_vec)). In the simplest model, each data point is a Bernoulli sample, but p is a function of time according to $$p(t) = e^{-t/T}$$. The goal is to estimate the time constant $$T$$.
It’s pretty easy to generate some fake data to practice with:
N_EXP = 200
time_const = 5
t_vec = np.linspace(0.1, 15, 20)
p_vals = np.exp(-1*t_vec/time_const)
data = np.random.binomial(n=1, p=p_vals, size=(N_EXP, len(p_vals)))
If I sum up the columns so the likelihood is Binomial, the sampling works just fine:
from theano.tensor import exp
with pm.Model() as binomial_model:
# Prior on time constant
BoundedNormal = pm.Bound(pm.Normal, lower=0.1)
tt = BoundedNormal('tt', mu=1., sd=10.)
pzero = pm.Deterministic('p', exp(-1*t_vec/tt))
obs = pm.Binomial('obs', n=N_EXP, p=pzero, observed=data.sum(axis=0))
trace = pm.sample(tune=1000, draws=5000)
It also works if I don’t sum and just use a Bernoulli likelihood with N_EXP independent trials:
with pm.Model() as bernoulli_model:
# Prior on time constant tt
BoundedNormal = pm.Bound(pm.Normal, lower=0.1)
tt = BoundedNormal('tt', mu=2., sd=10.)
# p decays in time
p = pm.Deterministic('p', exp(-1*t_vec/tt))
# Bernoulli likelihood with parameter p(t)
obs = pm.Bernoulli('obs', p=p, shape=len(t_vec), observed=data)
trace = pm.sample(tune=1000, draws=5000)
The more realistic situation is just a little more complicated, and this is where I can’t seem to get samples that make sense. Instead of measuring a bunch of 0 and 1 values directly, I have counts for each measurement which come from one of two Poisson distributions (parameterized by $$\mu_0$$ or $$\mu_1$$, where $$\mu_1 \neq \mu_0$$). The fake data now looks like this:
mu_0 = 2
mu_1 = 5
b = np.random.binomial(n=1, p=p_vals, size=(N_EXP, len(p_vals)))
data = (b*np.random.poisson(lam=mu_1, size=b.shape) + (1 - b)*np.random.poisson(lam=mu_0, size=b.shape))
My PyMC3 model with the Poisson likelihood:
with pm.Model() as poisson_model:
# Prior on time constant
BoundedNormal = pm.Bound(pm.Normal, lower=0.1)
tt = BoundedNormal('tt', mu=2., sd=10.)
# Form expected for decay of p(t)
pzero = exp(-1*t_vec/tt)
# Using p values draw Bernoulli samples
b = pm.Bernoulli('b', p=pzero, shape=len(t_vec))
# mu is either mu0 or mu1 depending on result of Bernoulli sample
muz = (b*mu_1 + (1 - b)*mu_0)
likelihood = pm.Poisson('like', mu=muz, shape=len(t_vec), observed=data)
trace = pm.sample(tune=1000, draws=5000)
Now I am not confident that this model is sampling properly. The distribution for tt is transformed and centers near the correct value, but if I look at the samples of b (the Bernoulli samples, parameterized by $$p(t)$$) the samples are always the same. The first few values are always one, and the rest are always zero. It seems like the sampler isn’t really moving around at all, at least for this variable. This is also reflected in the selection of $$\mu$$.
This behavior also shows up in the posterior predictive samples, which seem to always come from the $$\mu_1$$ Poisson distribution for the first few t values, then always from the $$\mu_0$$ distribution for the rest.
Any ideas why this isn’t working? I get the same behavior using pm.Categorical to select the Poisson parameters. The only model I have gotten reasonable looking PPC samples out of is a Poisson mixture model likelihood (with the columns of the data summed).
Thanks! I really appreciate any advice!