# Why does MCMC without observations have difficulty sampling?

**URL:** <https://discourse.pymc.io/t/why-does-mcmc-without-observations-have-difficulty-sampling/12809>\
**Category:** version agnostic\
**Tags:** prior, sampling\
**Created:** [August 28, 2023, 11:13am UTC](https://discourse.pymc.io/t/why-does-mcmc-without-observations-have-difficulty-sampling/12809 "2023-08-28T11:13:51Z")\
**Posts on this page:** 4\
**Page:** 1

<div class="post-metadata">

**Author:** ![Genesis](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/genesis/32/6017_2.png) [@Genesis](https://discourse.pymc.io/u/Genesis)\
**Post date:** [August 28, 2023, 11:13am UTC](https://discourse.pymc.io/t/why-does-mcmc-without-observations-have-difficulty-sampling/12809/1 "2023-08-28T11:13:51Z")

</div>

I’ve been trying to implement qa\_ver0 model by Michael Betancourt from [here](https://betanalpha.github.io/assets/case_studies/factor_modeling.html#52_Quality_Assurance) and found that removing the `observed` keyword from the likelihood specification massively slows down the sampling and causes significant sampling difficulties (divergences, very low effective sample size). Why is that?

I understand that with `observed` data, we are effectively sampling a prior. And I understand that the [recommended way to sample the prior](https://discourse.pymc.io/t/sampling-without-the-observed-keyword/10722) is via `sample_prior_predictive` (which uses [ancestral sampling instead of MCMC](https://discourse.pymc.io/t/pm-sample-how-can-the-posterior-be-defined-when-there-is-no-observed-data-to-calculate-the-likelihood/10930/5)). But I am not sure why MCMC sampling should be quite so ineffective. Yes it needs to convege, but it needs to converge with the posteriors too. And yes with a broad prior there is a large space for MCMC to explore - but still the sheer extent of the difficulties is surprising to me.

I think it’s a general question but for the sake of completeness here is my code. First with observations:

```auto
with pm.Model(coords={"observation_sample": range(N), "config": range(N_factor_configs)}) as qa_ver0:
    t_data = pm.ConstantData("obs_times", t, dims="observation_sample")
    y_data = pm.ConstantData("obs_passing_items", y, dims="observation_sample")
    
    baseline = pm.Normal("baseline", mu=-3, sigma=2)
    alpha = pm.Normal("alpha", mu=4, sigma=2)
    
    p = pm.Deterministic("p", pm.math.invlogit(alpha + baseline * t_data / 365.0), dims="observation_sample")
    
    passing_items = pm.Binomial("passing_items", p=p, n=N_samples, observed=y, dims="observation_sample")

with qa_ver0:
    trace_qa_ver0 = pm.sample()

```

No problems sampling:

```auto
Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 5 seconds.

```

Now if I remove `observed` data:

```auto
with pm.Model(coords={"observation_sample": range(N), "config": range(N_factor_configs)}) as qa_ver0_prior:
    t_data = pm.ConstantData("obs_times", t, dims="observation_sample")
    y_data = pm.ConstantData("obs_passing_items", y, dims="observation_sample")
    
    baseline = pm.Normal("baseline", mu=-3, sigma=2)
    alpha = pm.Normal("alpha", mu=4, sigma=2)
    
    p = pm.Deterministic("p", pm.math.invlogit(alpha + baseline * t_data / 365.0), dims="observation_sample")
    
    passing_items = pm.Binomial("passing_items", p=p, n=N_samples, dims="observation_sample")

with qa_ver0_prior:
    trace_qa_ver0_prior = pm.sample()

```

```auto
Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 920 seconds.
The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details
The effective sample size per chain is smaller than 100 for some parameters. A higher number is needed for reliable rhat and ess computation. See https://arxiv.org/abs/1903.08008 for details

```

To emphasise the key points here: 15 min to sample, high rhat, low effective sample size. In this specific example there are no divergences but previously I had also observed divergences

(For replication, the data (`y`, `t`, `N`, `N_samples` etc) are from [here](https://github.com/betanalpha/knitr_case_studies/blob/8c568575d8f4ca7700cf9b83faf73dc4a2089df2/factor_modeling/data/qa.data.R))

---

<div class="post-metadata">

**Author:** ![ricardoV94](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/ricardov94/32/5775_2.png) [@ricardoV94](https://discourse.pymc.io/u/ricardoV94)\
**Post date:** [August 28, 2023, 11:17am UTC](https://discourse.pymc.io/t/why-does-mcmc-without-observations-have-difficulty-sampling/12809/2 "2023-08-28T11:17:46Z")

</div>

I suggest plotting the prior draws and seeing how they look like

---

<div class="post-metadata">

**Author:** ![Genesis](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/genesis/32/6017_2.png) [@Genesis](https://discourse.pymc.io/u/Genesis)\
**Post date:** [August 28, 2023, 12:33pm UTC](https://discourse.pymc.io/t/why-does-mcmc-without-observations-have-difficulty-sampling/12809/3 "2023-08-28T12:33:40Z")

</div>

To be honest, I am not sure what I am looking for; hence the question, really…

Here is the plot of the prior:

```auto
az.plot_trace(trace_qa_ver0_prior.isel(observation_sample=0));

```

 ![prior_trace](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/1/19f21d69e0dc773be0d382ba35aeb2f5a4f8feb3.png)

Evidently, the `baseline` and `alpha` did not converge to Normal(-3,2) and Normal(4,2), respectively. But I guess we could have said that before too. There is hardly any mixing and chains are completely different - but again I am not sure how to link this to the the fact that we are sampling from a prior (since that doesn’t happen when I add `observed` and sample from the posterior)

There is definitely a problem with the fact that the `baseline` and `alpha` are then fed into a Binomial. If I remove the Binomial entirely then `pm.sample()` samples Normal(-3,2) and Normal(4,2) just fine:

```auto
with pm.Model(coords={"observation_sample": range(N), "config": range(N_factor_configs)}) as qa_ver0_minimal_prior:
    baseline = pm.Normal("baseline", mu=-3, sigma=2)
    alpha = pm.Normal("alpha", mu=4, sigma=2)
    
with qa_ver0_minimal_prior:
    trace_qa_ver0_minimal_prior = pm.sample()

az.plot_trace(trace_qa_ver0_minimal_prior.isel(observation_sample=0));

```

 ![prior_minimal_trace](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/b/b691556e2fe8ffde7f3dda7d409c44f32df7f673.png)

Conversely, if I fix `baseline` and `alpha` and sample only from the Binomial, again I have problems (although less extreme)

```auto
with pm.Model(coords={"observation_sample": range(N), "config": range(N_factor_configs)}) as qa_ver0_prior_minimal_binom:
    t_data = pm.ConstantData("obs_times", t, dims="observation_sample")
    y_data = pm.ConstantData("obs_passing_items", y, dims="observation_sample")
    
    alpha = -3
    baseline = 4
    
    p = pm.Deterministic("p", pm.math.invlogit(alpha + baseline * t_data / 365.0), dims="observation_sample")
    
    passing_items = pm.Binomial("passing_items", p=p, n=N_samples, dims="observation_sample")

with qa_ver0_prior_minimal_binom:
    trace_qa_ver0_minimal_binom_prior = pm.sample()

az.plot_trace(trace_qa_ver0_minimal_binom_prior.isel(observation_sample=0));

```

```auto
Sampling 4 chains for 1_000 tune and 1_000 draw iterations (4_000 + 4_000 draws total) took 550 seconds.
The rhat statistic is larger than 1.01 for some parameters. This indicates problems during sampling. See https://arxiv.org/abs/1903.08008 for details

```

 ![prior_minimal_binom_trace](https://canada1.discourse-cdn.com/flex036/uploads/pymc3/original/2X/e/ecb27d4d4ac027dd567713aa2738c8740df77327.png)

So that tells me that `pm.sample()` only recovers priors using NUTS but not Metroplis. Is that correct? If so, why? And why only for priors and not posteriors?

The second question that then arises - given that `pm.sample()` seems to deal fine with sampling normals via NUTS but doesn’t deal as well with the Binomial, why do the errors in the Binomial sampling “propagate back” to the Normals (and become exacerbated), even though there is no observed data - is that just a fundamental feature of MCMC? Is that one of the reasons why we use ancestral sampling for prior sampling - to “decouple” the variables in more complex models?

---

<div class="post-metadata">

**Author:** ![aseyboldt](https://yyz2.discourse-cdn.com/flex036/user_avatar/discourse.pymc.io/aseyboldt/32/5795_2.png) [@aseyboldt](https://discourse.pymc.io/u/aseyboldt)\
**Post date:** [August 28, 2023, 2:21pm UTC](https://discourse.pymc.io/t/why-does-mcmc-without-observations-have-difficulty-sampling/12809/4 "2023-08-28T14:21:07Z")

</div>

I guess another way of stating the question would be: If we forget all the structure that we know about the prior, and just treat is as a black-box with a likelihood (plus a bit of info about the domain of the prior), how hard is it to sample from that distribution?

I guess intuitively it might seem that the prior should be easier than the posterior, but I don’t think that really needs to be the case. It certainly can be, but it can just as easily be much harder.

Take for instance a simple model like this:

```auto
with pm.Model():
    sigma = pm.HalfNormal("group_sigma")
    group_effect = pm.Normal("group_effect", mu=0, sigma=sigma, dims="group")
    
    # Some likelihood that allows precise estimation of sigma:
    pm.Normal("y", mu=group_effect[group_idx], observed=some_large_dataset)

```

With the dataset, the prior might have a nice shape, and might be quite easy to sample. But without the data, this is a terrible funnel, that will not converge at all!

It looks like this is not really what’s happening in your example though. A big difference between sampling the prior and the posterior, is that in the first case pymc will also sample the dataset, but when you sample the posterior the dataset is assumed to be fixed. Since you are using a discrete likelihood, sampling the prior means sampling discrete parameters, where nuts doesn’t work. Sampling the posterior though doesn’t have that issue, and we can happily use gradient based methods.

But sampling distributions with discrete parameters is really hard. So I think it is not at all unexpected that this might fail horribly for even relatively small and simple (looking) problems.