I am trying to implement a Zero One Inflated Beta Regression model in PyMC3. Here’s what I have done so far,
with pm.Model(): # Priors for unknown model parameters theta1 = pm.Normal("theta1", mu=theta1_mu, sigma=theta1_sd) theta2 = pm.Normal("theta2", mu=theta2_mu, sigma=theta2_sd) theta3 = pm.Normal("theta3", mu=theta3_mu, sigma=theta3_sd) theta4 = pm.Normal("theta4", mu=theta4_mu, sigma=theta4_sd) theta5 = pm.Normal("theta5", mu=theta5_mu, sigma=theta5_sd) theta6 = pm.Normal("theta6", mu=theta6_mu, sigma=theta6_sd) precision = pm.Uniform("precision", 0, 100) # based on the original code # Expected value of outcome PI_0 = pm.math.invlogit(theta3_mu + theta4_mu * v) PI_1 = pm.math.invlogit(theta5_mu + theta6_mu * v) mu = pm.math.invprobit((1/theta2_mu) * pm.math.log(v/theta1_mu)) # Likelihood (sampling distribution) of observations y_zero = pm.Bernoulli("y_zero", p=PI_0, observed=y_zero_obs) y_one = pm.Bernoulli("y_one", p=PI_1, observed=y_one_obs) y_beta = pm.Beta("y_beta", alpha=mu * precision, beta=(1 - mu) * precision, observed=y_beta_obs) trace = pm.sample(1500, tune=1000, chains=3, cores=1) # tuned samples are being discarded
Since the probabilities PI_0, PI_1 are distributions, the model described here cannot be used as per my understanding.
The problem is that I have to mix these distributions (y_zero, y_one, y_beta) to have the following. Note that this is the likelihood.
However, with the code I have this is not happening. Secondly, I believe that the posteriors are not updating as per the likelihood, instead they are just recreating the priors. Could someone help me on this?
Thanks in advance!