I am trying to add a hierarchy on a constrained weight, but the centered parameterization is tripping me up. I think I have the non-centered version, but would prefer to test the centered version first as that is best practice.
# w number of weights, g number of groups
## non-centered
sigma_w = pm.HalfNormal('sigma_w', sigma=1, shape = (g))
z_w = pm.HalfNormal('z_w', sigma=1, shape = (w,g))
weight = pm.Deterministic('weight', z_w*sigma_w, shape=(w,g))
## centered?
sigma_w = pm.HalfNormal('sigma_w', sigma=1, shape = (g))
weight = pm.HalfNormal('weight', sigma=sigma_w, shape = (w,g))
Neither centered nor non-centered is “best practice”, there are just cases where you want to use one instead of the other. In general, non-centered is better when the sigma parameter admits quite small values (the groups are very similar), because this makes “funnels” that are hard to sample, even for NUTS (see here).
Otherwise, I don’t see anything wrong with your implementation. Did it somehow fail for you?
No. The hierarchy didn’t show much of a difference when added. I mainly want to make sure the centered version looks correct.
I would probably parametrize it in terms of a lognormal or something that has a concept of mean and dispersion. HalfNormal is limited. For instance, you can’t ever learn that most sigmas are tightly grouped around 3 (or any other non-zero value)
How would you do something like that? I agree on the dispersion aspect.
Very much like you would write a hierarchical normal prior but using lognormal or gamma for the mean and std
1 Like
This wouldn’t ensure the final output is strictly positive right? You would want to just take the exp of the output:
group_mu = pm.Normal('group_mu')
group_sigma = pm.HalfNormal('group_sigma', 1)
group_effect = pm.math.exp(pm.Normal('group_effect', mu=group_mu, sigma=group_sigma, dims=['group']))
Could also take abs to get a “generalized folded normal”
That’s why I mentioned lognormal, which has that implicit exp baked in
3 Likes