Generating population level samples for heirarchical parameters

Questions
1

In one of my models, I estimate parameters for several subjects through the following way:

    log_param_mean = pm.Bound(pm.Normal, lower=0)("log_param_mean", tt.log(3.5), 1)
    rfx = pm.Normal("rfx", 0, 1, shape=len(np.unique(subject_ids)))
    sigma_param = pm.HalfCauchy("sigma_param", 1)

    params = pm.Deterministic("params", tt.exp(log_param + rfx[np.unique(subject_ids)] * sigma_param))

If I were to sample the posterior predictive distribution with PyMC3, then I would get posterior samples for each of my subjects. What I’m really interested in is population level samples. So right now, this parameter has a log-normal distribution. How can I use PyMC3 to sample that log-normal distribution and get population level samples for my parameter (either when using pm.sample or pm.sample_posterior_predictive)?

Could I do something like

 rng = tt.shared_randomstreams.RandomStreams()
sampled_param = pm.math.exp(log_param_mean + rng.normal()*sigma_param)
pop_param = pm.Deterministic('pop_param', sampled_param)

I hope that is clear. If not, I can gladly add context.

2

Hi Demetri!
I’m not sure I completely understand your issue, but yeah your proposal looks good to me. Then you can use var_names inside pm.sample_posterior_predictive to sample pop_param too – instead of only the observed variable.
Hope this helps :vulcan_salute:

3

Hey Alex,

I think what I’m trying to do is something similar to Stan’s generated quantites block. I want to use my samples to make new quantities. Looking for the safest way to do that.

4

Mmmh I’d say pm.sample_posterior_predictive would be the way to do that, or just pushing your parameters through the model, but maybe I don’t understand what you’re trying to do :thinking:
Perhaps @ahartikainen would have valuable inputs here?

5

Hi,

One way to do this is using e.g. ArviZ and calculating your values after the posterior has been sampled.

I don’t know if pm.Deterministic could work in this case.

2 Likes
6

“calculating your values after the posterior has been sampled”: That’s what I meant by “pushing your parameters through the model” – thanks Ari for putting well-defined words on my confused thoughts :laughing:

And I’d say pm.Deterministic would work then, as it’s just a wrapper to compute determined combinations of posterior parameter values.