Adding levels to a model with Gamma & Exponential priors

Misc.
1

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)

2

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

3

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

Think I got it :slight_smile:

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