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)
```