I am working on a modeling problem where the regression is perfectly suited for a hierarchical structure. However, after inspecting the data, it seems that the response variable has different mean and standard deviation for each group. This might be a straightforward problem, but I wondering how to extend the multilevel structure to the noise in the likelihood. e.g., if we have some x-data and some y-data with four groups, we might model like likelihood as:

```
intercept = pm.Normal(..., dims=groups)
slope = pm.Normal(..., dims=groups)
mu = intercept[group] + slope[group]*x-data
err = pm.HalfNormal(...)
likelihood = pm.Normal(..., mu = mu, sigma = err, dims="obs_id")
```

This works fine in general, but it assumes a single noise param for all groups. Is it as simple as changing the noise to be something like:

```
err = pm.__SOME_DISTRIBUTION__('err', .... , dims=groups)
```

I wrote a full working example with simulated data. This question is based off of a problem that I am actually working on, but for simplicity in understanding, I am working with the simulated set. I was able to get the above approach to work, but (a) I am not sure if this is a correct approach, and (b) how to appropriately model the noise when I need prior and hyperprior.

Here is the notebook on colab. I thought it was clearer than trying to fit all the code/figures into one long post.

Thanks!