Hello everyone!

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!