Many thanks @jhrcook and @OriolAbril for the answers! Only in the last days I managed to find the time to get back to this, and the solutions you suggested indeed allowed me to calculate HDIs.
I’m now facing another problem though, my predictions using the sample_posterior_predictive function seem to be way off the expected target (by an order of magnitude approximately). On the other hand, if instead of using the sample_posterior_predictive function, I take the parameters directly from the trace, the predictions seem to be pretty accurate.
More concretely, if I want to calculate predictions for the model that I described in the original post, the following approach:
with partial_pooling:
pm.set_data({"daypart":daypart, “temperature”: temp_ts})
posterior_predictive = pm.sample_posterior_predictive(partial_pooling_trace, var_names=[“a_daypart”, “bt_daypart”])
a_daypart_means_1 = posterior_predictive[‘a_daypart’].mean(0)
bt_daypart_means_1 = posterior_predictive[‘bt_daypart’].mean(0)
predictions_1 = a_daypart_means_1[np.ix_(dayparts)] + bt_daypart_means_1[np.ix_(dayparts)] * temp_ts
Produces predictions that are way off target, while:
a_daypart_means_2 = np.mean(partial_pooling_trace[‘a_daypart’], axis =0)
bt_daypart_means_2 = np.mean(partial_pooling_trace['bt_daypart’], axis =0)
predictions_2 = a_daypart_means_2[np.ix_(dayparts)] + bt_daypart_means_2[np.ix_(dayparts)] * temp_ts
produces pretty accurate predictions. I’m quite sure that the first way is the correct way to proceed, even though I’m not really aware of the difference between the two approaches.
Could someone shed some light on that? Is the second approach causing severe overfitting? Also if anyone could suggest some ideas to improve the predictions obtained with the first approach, that would be amazing.
Thanks a lot again for the great support provided on Pymc discourse!