I’m trying to extend the hierarchical regression model in this tutorial for a multi-dimensional dataset (dataset with multiple features).
I was trying to change the original model as shown below.
with pm.Model() as hierarchical_model:
# Hyperpriors for group nodes
mu_a = pm.Normal('mu_a', mu=0., sd=100**2)
sigma_a = pm.HalfCauchy('sigma_a', 5)
mu_b = pm.Normal('mu_b', mu=0., sd=100**2)
sigma_b = pm.HalfCauchy('sigma_b', 5)
# Intercept for each county, distributed around group mean mu_a
# Above we just set mu and sd to a fixed value while here we
# plug in a common group distribution for all a and b (which are
# vectors of length n_counties).
a = pm.Normal('a', mu=mu_a, sd=sigma_a, shape=n_counties)
# Intercept for each county, distributed around group mean mu_a
b = pm.Normal('b', mu=mu_b, sd=sigma_b, shape=(n_counties, 2))
# Model error
eps = pm.HalfCauchy('eps', 5)
radon_est = a[county_idx] + tt.dot(b[county_idx], data[["floor", "wave"]].values)
# Data likelihood
radon_like = pm.Normal('radon_like', mu=radon_est, sd=eps, observed=data.log_radon)
with hierarchical_model:
hierarchical_trace = pm.sample(draws=2000, n_init=1000)
However, I get an error when trying to perform the dot product between the random variable b and the predictors (X). Here, instead of single predictor “floor” I have selected the feature “wave” in order to test the model with multiple features.
Can someone please help me to figure out the correct way to implement this model
Thanks !!