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
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.
I’m not sure I completely understand your issue, but yeah your proposal looks good to me. Then you can use
pm.sample_posterior_predictive to sample
pop_param too – instead of only the observed variable.
Hope this helps
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.
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
Perhaps @ahartikainen would have valuable inputs here?
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.
“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
And I’d say
pm.Deterministic would work then, as it’s just a wrapper to compute determined combinations of posterior parameter values.