Thanks again for your answers @OriolAbril @jhrcook. I carefully read the posts and the notebook mentioned. So, from what I understand, sampling variables of the posterior to obtain predictions does not make sense. Instead, the following methods should provide the correct predictions:
1- sampling the posterior distribution using pm.sample_posterior_predictive(partial_pooling_trace)["y"].
2 - manually extracting the individual parameter means from the trace, like I did in the second approach I previously posted, and then compute each prediction using predictions_2 = a_daypart_means_2[np.ix_(dayparts)] + bt_daypart_means_2[np.ix_(dayparts)] * temp_ts
3- another option might be manually modifying the data to include all the possible independent variable combinations once, as suggested in Oriol’s notebook, but since one of the independent variables I have is a long timeseries of continuous data, I think this option might not work well in my case (seems more oriented to categorical data).
I understand that the first option should be the most ‘correct’? Or is the second one a valid alternative as well?
Thanks again for the enlightening comments and links!