I ended up doing exactly this as in my use-case (a rather large model) it ended up being faster sampling everything together (3h30) than doing it sequentially using the empirical distribution (5h30) discussed in these other threads and this pymc-experimental issue. However that might be due to implementation inefficiencies of the proposed empirical distribution, or something else. (I used pymc 5.9 with python 3.11 because the model breaks with subsequent pymc versions for some reason.) I keep thinking a proper (optimized) tool for this would be great to have in a toolbox. It would certainly be great for teaching Bayesian statistics !
2 Likes