# Bernoulli samples all the same?

**URL:** <https://discourse.pymc.io/t/bernoulli-samples-all-the-same/1297>\
**Category:** Questions\
**Created:** [May 29, 2018, 6:21pm UTC](https://discourse.pymc.io/t/bernoulli-samples-all-the-same/1297 "2018-05-29T18:21:28Z")\
**Posts on this page:** 2\
**Page:** 1

<div class="post-metadata">

**Author:** ![bod](https://avatars.discourse-cdn.com/v4/letter/b/ed655f/32.png) [@bod](https://discourse.pymc.io/u/bod)\
**Post date:** [May 29, 2018, 6:21pm UTC](https://discourse.pymc.io/t/bernoulli-samples-all-the-same/1297/1 "2018-05-29T18:21:28Z")

</div>

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:

```python
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:

```python
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:

```python
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:

```python
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:

```python
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$$.

![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/1X/0236cb0bf10307b5856fc381d9160b14a0c404b8.png)

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.

![image](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/1X/ffa3d3d663f0d476665a231750f7c302a293f033.png)

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!

---

<div class="post-metadata">

**Author:** ![junpenglao](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/junpenglao/32/8_2.png) [@junpenglao](https://discourse.pymc.io/u/junpenglao)\
**Post date:** [May 29, 2018, 7:10pm UTC](https://discourse.pymc.io/t/bernoulli-samples-all-the-same/1297/2 "2018-05-29T19:10:27Z")

</div>

Metropolis usually perform poorly in high dimension - and in this case, the proposal value during sampling the Bernoulli value got rejected - which gives a stuck MCMC chain.

For example, here is a similar model with the same problem, and I did some diagnostic:

> [@Indexing using a free random variable](https://discourse.pymc.io/t/indexing-using-a-free-random-variable/939/15):
>
> I had another look at the model, with the following code: with pm.Model() as basic\_model: # Priors for unknown model parameters b1 = pm.Uniform('b1', lower=0.3, upper=0.5, testval=0.45) ncomp\_aCH = pm.Categorical('ncomp\_aCH', p=np.ones(n)/n) ncomp\_aCOH = pm.Categorical('ncomp\_aCOH', p=np.ones(n)/n) aCH=aCH\_[0, ncomp\_aCH] aCOH=aCOH\_[0, ncomp\_aCOH] out= b1\*aCH+aCOH # Likelihood (sampling distribution) of observations Y\_obs = pm.Normal('Y\_obs', mu=out, tau=si…

The best solution is to rewrite your model into a mixture model.