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