Hi everyone,
I’d like to better understand pymc3 and hierarchical model. I am currently stuck on how to properly estimate the following problem, which - in my mind - should be solvable via a hierarchical structure (but maybe it isn’t?). Looking for help but not necessary an answer, a pointer to some tutorial that I may have missed would be helpful.
TLDR
Key questions in the following wall of text
1) How do I append priors for the different levels?
2) How do I append data for the posterior sampling for the different levels? (and can I do that)
Below it is broken down in problem setup, (pseudo) code, key-questions on coding practices/hierarchical models
# Problem Setup
Summary
I have sales information for a product at a National, State, County, etc… level.
The model, true for each level, I am trying to fit is the following:
log(Y_l,i,j) = a_l,i + b_0_l,ixPrice_l,i,j + b_1_l,ixDiscount_l,i,j + b_2_l,ixWeekOfMonth_l,i,j
Where:
l = level
i = group within that level (e.g. California)
j = daily observation for that group
The b_0, b_1 variables should be hierarchical, meanwhile the others shouldn’t. Reason for this lies in the fact that the consumer response to price and discount needs to be similar and “consistent” to each other, meanwhile the other variables are significantly more regional.
To visualize, my data looks something like the following:
| Week Of Month | Level | Units | Const | Edp | Discount |
|---|---|---|---|---|---|
| 1 | USA | 10000 | 1 | 0.5 | 0 |
| 1 | Alabama | 100 | 1 | 0.3 | 0.01 |
| 1 | Arizona | 200 | 1 | 0.4 | 0.02 |
| 1 | Cook county | 10 | 1 | 0.22 | 0.04 |
| 1 | Elmore County | 10 | 1 | 0.321 | 0 |
Given what said before, I would like to set the following priors:
for each a_l,i = Normal(mean(log(Y_l,i,k)), std(log(Y_l,i,j)))
for each b_2_l,i = InverseGamma ()
For level 0:
b_0_0 = Normal
For level_1:b_0_1,i = Normal(b_0_0)
For level_2:b_0_2,k = Normal(b_0_1,i)
Where k in this case represent the county belonging to state i
# Pseudo Code
Summary
Looking at the documentation, I have written the following “pseudo” code (I’ll highlight why i call it pseudo-code):
with pm.Model() as varying_intercept_slope:
# Nation Priors
constant_nation = pm.Lognormal('constant_nation', p.mean(log_gross_units_nation), np.std(log_gross_units_nation)))
edp_nation = pm.Normal('mu_a', mu=0., sigma=2)
sigma_edp_nation = pm.HalfCauchy('sigma_a', 5)
discount_nation = pm.Normal('mu_b', mu=0., sigma=2)
sigma_discount_nation = pm.HalfCauchy('sigma_b', 5)
week_nation = pm.Normal('week_nation', mu = 0, sigma = 1)
# State Priors
# Random constant
constant_states = []
sigma_edp_state = []
sigma_discount_nation = []
for i in range(0,50):
constant_states.append(pm.Lognormal(np.mean(log_gross_units_states[i]), np.std(log_gross_units_state[i])))
sigma_edp_state.append( pm.HalfCauchy('sigma_a', 5) )
sigma_discount_nation.append( pm.HalfCauchy('sigma_b', 5) )
# Random edp states
edp_states = pm.Normal('edp_states', mu = edp_nation, sigma=sigma_edp_nation, shape=states)
# Random discount states
discount_states = pm.Normal('discount_states', mu=discount_nation, sigma=sigma_discount_nation, shape=states)
# Random week coeffs
week_states = pm.Normal('week_states', mu = 0, sigma = 1, shape = states )
# Counties Priors
# Random constant
constant_counties = []
edp_counties = []
discount_counties = []
week_counties = []
week_states.append(('week_states', mu = 0, sigma = 1, shape = counties ))
for i in range(0,len(counties)):
constant_counties.append(pm.Lognormal(np.mean(log_gross_units_counties[i]), np.std(log_gross_units_counties[i])))
for i in range(0,50):
# Random edp states
edp_counties.append( pm.Normal('edp_states', mu = edp_states[i], sigma=sigma_states[i], shape=states))
# Random discount states
discount_states = pm.Normal('discount_states', mu=discount_states[i], sigma=sigma_discount_states[i], shape=countie)
# Model error
sigma_y = pm.Uniform('sigma_y', lower=0, upper=100)
# Expected value
y_hat_nation = constant_nation + edp_nation * edp + discount_nation * discount + week_nation * week
y_hat_states = constant_state + edp_state[state] * edp + discount_state[state] * discount + week_county[county] * week
y_hat_counties = constant_counties + edp_county[county] * edp + discount_county[county] * discount + week_county[county] * week
y = [y_hat_nation, y_hat_states, y_hat_counties]
# Data likelihood
y_like = pm.Normal('y_like', mu = y, sigma=sigma_y, observed=gross_units)
Key code questions
In a situation like this, where I am recursively setting priors, how do I append the priors to each other?
for i in range(0,50):
# Random edp states
edp_counties.append( pm.Normal('edp_states', mu = edp_states[i], sigma=sigma_states[i], shape=states))
Since I am trying to forecast the gross units at all three levels, how do I write down the equation/how do I append the dependent variables to each other?
# Expected value
y_hat_nation = constant_nation + edp_nation * edp + discount_nation * discount + week_nation * week
y_hat_states = constant_state + edp_state[state] * edp + discount_state[state] * discount + week_county[county] * week
y_hat_counties = constant_counties + edp_county[county] * edp + discount_county[county] * discount + week_county[county] * week
y = [y_hat_nation, y_hat_states, y_hat_counties]
# Data likelihood
y_like = pm.Normal('y_like', mu = y, sigma=sigma_y, observed=gross_units)
Is it possible to have a mixed hierarchical model, were some of the coefficients are not hierarchical? i.e. how do I set my constant to be a free parameter from a distribution
Thanks for the help
