Hello PyMC3 community,
I’m coming to you with another question about hierarchical modeling for regression. My problem is that sampling is taking way too long. When I try to sample 2000 draw on 4 chains with 500 tuning steps, it takes nearly 5 hours to run, which will make my future plan of cross-validation take an unreasonable amount of time.
I believe the reason that my sampling is taking so long is related to the hierarchal structure I placed on my intercept term. My idea was that my intercept term should be informed by county level, state level and country level observations; something along the lines of:
a_county = a_state_mean + county_level_noise
a_state = a_country_mean + state_level_noise
a_country = global_mean + noise
where the idea is that if I get an observation from a new county, then I should at least have the state mean as a hopefully good enough starting point for making a prediction. The code for my model is as follows:
with pm.Model(coords=coords) as work_in_progress: # Hyper Priors: Country sigma_a_country = pm.HalfNormal("sigma_a_country", 3.) mu_a_country = pm.Normal("mu_a_country", mu=9., sigma=2.5) # Country coefficients a_country_offset = pm.Normal("a_country_offset", mu=0., sigma=1., dims= "country") a_country = pm.Deterministic("a_country", mu_a_country + sigma_a_country*a_country_offset) # Hyper Priors: State sigma_a_state = pm.HalfNormal("sigma_a_state", 3.) # State coefficients a_state_offset = pm.Normal("a_state_offset", mu=0., sigma=1., dims= "state") a_state = pm.Deterministic("a_state", a_country[country] + sigma_a_state*a_state_offset[state]) # Hyper Priors: County sigma_a_county = pm.HalfNormal("sigma_a_county", 3.) # County coefficients a_county_offset = pm.Normal("a_county_offset", mu=0., sigma=1., dims= "county") a_county = pm.Deterministic("a_county", a_state[state] + sigma_a_county*a_county_offset[county]) # Model error epsilon = pm.HalfCauchy("eps", 5.0) target_estimate = a_county[county] # Data likelihood target_like = pm.TruncatedNormal("target_likelihood", mu=yield_estimate, sigma=epsilon, observed = target.values, dims="obs_id", lower = 0.)
When I visualize this model
pm.model_to_graphviz(work_in_progress) I think the issue is that both
a_county are at the observation level and I’m not sure how to fix that.
I’m new to PyMC3/Bayesian modeling so maybe my thought process is wrong but in my mind, they should both be like
a_country. Meaning they should both be contained in their own box of dimension 35 and 187 respectively and the only thing that should be at the observation level is the likelihood?
Any feedback would be greatly appreciated. The closest similar question I could find was here but I’m not sure how this individual fixed their intercept terms. Any feedback would be greatly appreciated.