I need help thinking through my approach to Gibbs Sampling of many parameters and I'd like to know if there is literature on this topic:
I have a dataset with 3 dimensions: state, year, code. With the terms and interaction terms, there are $$50 + 5 + 5 + 50*5 + 5*5 + 50*5*5 = 1,585$$ parameters that I need to estimate. I am using Gibbs Sampling (JAGS via R) to fit a model to the data. I am running 3 chains (each on a separate processor). With 1,585 parameters to estimate, JAGS takes about 12 hours to run.
Now, I want to swap the code dimension for another dimension: type. (Note that type and code are separate fields and type is not a refinement of code.) I am also removing the state dimension. This results in $$5 + 750 + 5*750 = 4,505$$ parameters to estimate.
From testing, it appears that each added parameter takes exponentially more time to fit. So this model cannot complete within a week.
I have access to several servers (about 100 cores in total), so I want to take advantage of that. The idea: instead of doing 1 run of Gibbs Sampling with all parameters, do $n$ runs of Gibbs Sampling in parallel on subsets of parameters.
Process:
- Perform stratified sampling to select a subset of values from
type(sample 100 values, instead of all 750 at once). - All unselected data is aggregated to create an 'all other' bucket. I am doing this so that the population means across all runs are similar.
- Run the model 200 times (each time estimating 605 parameters which were selected prior to the run).
- At the end, average the estimates that were generated in the runs of Gibbs Sampling (so that there is only 1 estimate per parameter).
This makes me think of Random Forests in modeling decision trees (where you sample from the dataset, run the model on a subset of data and average the results), but that is in the context of empirical models. Does similar logic work for Bayesian hierarchical models? Do you see a flaw in this approach? And is there literature on this topic? Thanks!
