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”