Thank you very much for the answer Thomas. I was indeed able to calculate my predictions by using the pm.sample_posterior_predictive() function as described in the notebook you linked.
I’m now trying to calculate HDIs for my predictions using the az.hdi() function, but I’ve been encountering different issues. I would like to be able to estimate an HDI interval for my predictions, and my idea is that in order to obtain that I should calculate HDIs for each of the posterior distributions of the parameters in my model, and then combine them as specified by the model. But I’m not sure if this is the correct way, and I’m also not sure of how to obtain said HDIs, since when I try to run for example
az.hdi(posterior_predictive['a_daypart'])
I get the error:
UserWarning: More chains (4000) than draws (6). Passed array should have shape (chains, draws, *shape)
which I think is related to the fact that my posterior predictive is really just 1 chain with 4000 draws, while 6 is the amount of different ‘dayparts’, but arviz doesn’t get that.
Overall I’m just really confused about how the arviz data structure works and how to calculate HDIs in general and I’d be really grateful if someone could shed some light or point to any resources that can help me!