How do I put a prior on a deterministic/intermediate variable?

I hope this question makes sense. I’m trying to reproduce something from the paper Bayesian Analysis of Errors-in-Variables Regression Models (Dellaportas 1995). They model the growth rates of fish from tanks with different levels of dissolved oxygen. They have a slope \beta_i for each tank i, which is defined as the deterministic function:

\beta_i = \frac{O_{ai} - \phi_1}{(\frac{O_{ai}}{\phi_2}) - \phi_3}

where O_{ai} is a latent variable for the true oxygen level in tank i, \phi_1 and \phi_2 are parameters relating oxygen to the growth rate that are given informative priors, and \phi_3 is a parameter with an improper flat prior over the whole real line.

However, they also put priors on the \beta_i terms themselves. I think this has the effect of indirectly constraining \phi_3 so it’s not really an improper prior, but I’m not sure.

In any event, I’m not sure how to represent such a thing in pymc3. The beta terms have to be written as a function of the other variables, but they also need a prior. The only thing I can think to do is to declare beta as a deterministic result of the other variables, and then “bolt on” a prior by adding a pm.Potential term which penalizes the likelihood in the right way. Thoughts?

Here’s the full context from the paper:

Here’s some nonfunctioning code to give you an idea of where I’m at right now:

with pm.Model() as the_model:
    alpha = pm.Normal("alpha", mu=150, sd=50) # common intercept
    
    beta_means = np.array([9, 9, 9, 6.9, 4.6, 3.3])
    betas = pm.Normal("betas", mu=beta_means, sd=np.ones(6)) # tank prior slopes
    
    phi_1 = pm.Lognormal("phi_1", mu=0.28, sd=0.59) # critical oxygen value
    phi_2 = pm.Lognormal("phi_2", mu=1.73, sd=0.63) # saturation oxygen value
    phi_3 = pm.Flat("phi_3") # ?
    
    tau_w = pm.Gamma("tau_w", alpha=3.5, beta=2) # prior on e_w precision
    tau_b = pm.Gamma("tau_b", alpha=1, beta=1) # prior on e_b precision
    tau_o = pm.Gamma("tau_o", alpha=2, beta=1) # prior on e_o precision

    true_oxygen = pm.Flat("true_oxygen", shape=6)
    observed_oxygen = pm.Normal("observed_oxygen", mu=true_oxygen, tau=tau_o,
                                observed=oxygen_values)

    ############################################################################
    ### something with betas as a function of the phis and true_oxygen here? ###
    ############################################################################

    weight_means = alpha + betas*np.arange(1,9)[:, None]
    weights_observed = pm.Normal("weights_observed", mu=weight_means,
                                 tau=tau_w, observed=weight_values)

That’s how I would do as well.

1 Like

I haven’t looked at all the context of the problem, but I have been in a similar situation with pymc2. I wanted to apply a prior-distribution to a deterministic variable.

Fortunately, I had symmetries on my side. I made up a dummy data-value with the central value of my prior distribution. Then I told pymc2 that the dummy-data varied according to the deterministic variable. And in the wash, I was able to show that this was the same as applying a prior I wanted to use. That may or may not work for you.