Constraint or prior dependent on multiple parameters

Questions
1

I have a question about estimating a linear model with constraints in pymc3:

y = \beta_0 + \beta_1p + \beta_2t +\epsilon

Estimating this model without any other conditions is relatively straight forward, however I have prior knowledge on the ratio of two parameters:

\frac{\beta_1}{\beta_0} \sim N(\mu,\sigma^2)

Is the most direct way to implement a constraint? I do not know the individual distributions of either \beta_0 or \beta_1 so I cannot separate them.

Not sure if this is possible, or an entirely valid modelling approach so happy to be told its not possible if that is the case.

2

Usually, we express hard constraint using potential, as it effectively cut out some part of the parameter space.
In this case, you can do it in a similar way:

with pm.Model() as model:
    inv_b0 = pm.Flat('inv_b0')
    b1 = pm.Flat('b1')
    b0b1 = pm.Deterministic('b1/b0', b1*inv_b0)
    pm.Potential('constraint', pm.Normal.dist(5., 1.).logp(b0b1))
    trace = pm.sample()
pm.traceplot(trace, compact=True);

However, this model is unidentifiable as is, since you can have b0*, b1* = k*b0, k*b1 and still have b1* / b0* == b1 / b0. You would need stronger prior on b0 and b1 to make it work nicely.

Also, note that the model above is expressing prior knowledge of the ratio, it’s NOT the same as constraint (i.e., the condition that the posterior of the parameter also satisfy). To express this as constraint is a bit more difficult, I would go with something like:

mu, sigma = 5., 5.
rand_stm = tt.shared_randomstreams.RandomStreams()

with pm.Model() as model:
    b0 = pm.Normal('b0')
    b1_div_b0 = pm.Deterministic('b1/b0', rand_stm.normal(avg=mu, std=sigma))
    b1 = pm.Deterministic('b1', b1_div_b0*b0)
#     pm.Potential('b1_prior', pm.Normal.dist(5., 2.).logp(b1))
    trace = pm.sample()
pm.traceplot(trace, compact=True);

Note that now you actually have b1 / b0 ~ Normal(mu, sigma) even for the posterior.

2 Likes
3

Ok so in your first example is that prior completely useless, even if I have observational data to estimate the parameters?

Also is there a reason why you should define \frac{1}{\beta_0} as opposed to directly defining the parameters, and then b0b1 = pm.Deterministic('b1/b0', b1/b0)

I was thinking something like the below model:

npt = 100
x = np.linspace(0,200, npt)
y = 10 + 5 * x  + np.random.normal(0, 5, npt)

if __name__ == '__main__':
    with pm.Model() as model:
        sigma = pm.HalfNormal('sigma', sigma=20)
        b0 = pm.Flat('b0')
        b1 = pm.Flat('b1')
        b0b1 = pm.Deterministic('b1/b0', b1/b0)
        pm.Potential('constraint', pm.Normal.dist(1, 1).logp(b0b1))
        likelihood = pm.Normal('y', mu=b0 + b1 * x + b2 * z,
                            sigma=sigma, observed=y)
        trace = pm.sample(2000, tune=500, chains=2)
    pm.traceplot(trace, compact=True);

Which I believe you are saying is equivalent to setting a prior on the ratio of the parameters?

4

b1/b0 likely will give you invalid result as b0 could be 0.
The prior still plays a role here, and will be affected by the observation.

5

Thanks! seems to be working as intended, so marked your original response as the solution.

1 Like