I’m trying to fit a hierarchical model using data from 30 features and 5 groups (code below). Since I do not want to define each of the 30 features’ distributions, I’m using the
shape argument to speed up things. This will require to create a matrix of size
(n_groups, n_features) to store the coefficients of each feature per group. I’ve two questions:
- When I run this, I get an error:
Exception: ("Compilation failed (return status=1). I guess I’m doing something wrong when defining the model, but I cannot figure out what’s the mistake.
- My understanding is that the code below defines a model that assumes that the coefficients are independent of each other. If I wanted to consider dependence among coefficients, what would need to be modified?
Any reference / guidance for these questions is highly appreciated.
# X is a matrix of size (N, 30) # y is a vector of length N # group_id is a vector of length N with n_groups unique values n_features = 30 n_groups = 5 with pm.Model() as hierarchical_model: # Hyperpriors for group nodes mu_a = pm.Normal('mu_a', mu=0., sigma=10) sigma_a = pm.HalfNormal('sigma_a', sigma=20) mu_b = pm.Normal('mu_b', mu=0., sigma=10, shape=n_features) sigma_b = pm.HalfNormal('sigma_b', sigma=20, shape=n_features) a_offset = pm.Normal('a_offset', mu=0, sd=1, shape=n_groups) a = pm.Deterministic("a", mu_a + a_offset * sigma_a) b_offset = pm.Normal('b_offset', mu=0, sd=1, shape=(n_groups, n_features)) b = pm.Deterministic("b", mu_b + b_offset * sigma_b) # Model Error noise = pm.Gamma('noise', alpha=2, beta=1) estimate = a[group_id] + X @ b[group_id] y_observed = pm.Normal('y_observed', mu=estimate, sigma=noise, observed=y)