# Adding levels to a model with Gamma & Exponential priors

I have a hierarchical model that’s currently setup like this:

def model(x, y, level1_id, level2_id):
level0_intercept = numpyro.sample('level0_intercept', dist.Normal(0, 1))
level0_intercept_scale = numpyro.sample('level0_intercept_scale', dist.Exponential(1))

with numpyro.plate('level1', len(set(level1_id)), dim = -1):
level1_intercept = numpyro.sample('level1_intercept', dist.Normal(level0_intercept, level0_intercept_scale))
level1_intercept_scale = numpyro.sample('level1_intercept_scale', dist.Exponential(1))

with numpyro.plate('level2', len(set(level2_id)), dim = -2):
level2_intercept = numpyro.sample('level2_intercept', dist.Normal(level1_intercept, level1_intercept_scale))

intercept = level2_intercept[level2_id, level1_id]
slope = numpyro.sample('slope', dist.Gamma(2, 1))
curvature = numpyro.sample('curvature', dist.Exponential(1))

loc = intercept - slope*x - curvature*x**2
scale = numpyro.sample('error', dist.HalfNormal(1))

with numpyro.plate('observation', len(x)):
numpyro.sample('y_est', dist.Normal(loc, scale), obs = y))


What I want to do is to include slope & curvature in the hierarchy, just like the intercept.

But how do I do this for a gamma & exponential distribution?
For Gamma I guess I might be able to model concentration & rate in the hierarchy, but previous attempts has resolved in invalid values (<=0)

Similar to intercept_scale, you should use a positive support prior (e.g. HalfNormal, Exponential, LogNormal,…) for positive parameters.

Yes, but how do I connect those to the next level?

With a Normal() I add a scale distribution that’s unique for that level, e.g. level0_intercept_scale above, and use that as priors on the next level.

Do I just chain them?

level0_curvature = dist.Exponential(1)
level1_curvature = dist.Exponential(level0_curvature)
level2_curvature = dist.Exponential(level1_curvature)


Think I got it

For Gamma():
mean = \frac{concentration}{rate}
concentration = mean \cdot rate

That means I can take the output of the previous Gamma(), sample rate and then calculate concentration on the next level.

As for Exponential() & HalfNormal(), I took the inverse of the previous output as the input for the next level.

def model(x, y, level1_id, level2_id):
level0_intercept = numpyro.sample('level0_intercept', dist.Normal(0, 1))
level0_intercept_scale = numpyro.sample('level0_intercept_scale', dist.Exponential(1))

level0_slope = numpyro.sample('level0_slope', dist.Gamma(2, 1))
level0_slope_rate = numpyro.sample('level0_slope_rate', dist.Exponential(1))

level0_curvature = numpyro.sample('level0_curvature', dist.Exponential(1))

with numpyro.plate('level1', len(set(level1_id)), dim = -1):
level1_intercept = numpyro.sample('level1_intercept', dist.Normal(level0_intercept, level0_intercept_scale))
level1_intercept_scale = numpyro.sample('level1_intercept_scale', dist.Exponential(1))

level1_slope = numpyro.sample('level1_slope', dist.Gamma(level0_slope*level0_slope_rate, level0_slope_rate ))
level1_slope_rate = numpyro.sample('level1_slope_rate', dist.Exponential(1))

level1_curvature = numpyro.sample('level1_curvature', dist.Exponential(1/level0_curvature))

with numpyro.plate('level2', len(set(level2_id)), dim = -2):
level2_intercept = numpyro.sample('level2_intercept', dist.Normal(level1_intercept, level1_intercept_scale))
level2_slope = numpyro.sample('level2_slope', dist.Gamma(level1_slope*level1_slope_rate, level1_slope_rate))
level2_curvature = numpyro.sample('level2_curvature', dist.Exponential(1/level1_curvature))

intercept = level2_intercept[level2_id, level1_id]
slope = level2_slope[level2_id, level1_id]
curvature = level2_curvature[level2_id, level1_id]

loc = intercept - slope*x - curvature*x**2
scale = numpyro.sample('error', dist.HalfNormal(1))

with numpyro.plate('observation', len(x)):
numpyro.sample('y_est', dist.Normal(loc, scale), obs = y))

1 Like