I have a model (simplified model below) that has **multiple observed variables** that depend on each other. After model training, I’d like to predict by sampling from the posterior.

It seems like this version does not allow me to successfully predict both target variables.

Generate some example data:

```
seed = 42
np.random.seed(seed)
def get_data(n = 100):
mu = 1000
X_obs = stats.poisson.rvs(mu=mu, size=n).astype(float)
r_obs = stats.beta.rvs(a=100, b=50, size=n)
br_obs = stats.beta.rvs(a=90, b=10, size=n)
w_obs = [0.25, 0.75]
y_obs = X_obs * r_obs * br_obs
return y_obs, r_obs, br_obs, w_obs, X_obs
n = 1000
y_obs, r_obs, br_obs, w_obs, X_obs = get_data(n)
R_obs = r_obs * X_obs
```

```
X_shared = tt.shared(X_obs)
```

```
with pm.Model() as model_not_predictable:
sd_alpha = pm.HalfNormal("sd_a", sd=1, shape=2)
sd_beta = pm.HalfNormal("sd_b", sd=1, shape=2)
alpha = pm.HalfNormal("alpha", sd=sd_alpha, shape=2)
beta = pm.HalfNormal("beta", sd=sd_beta, shape=2)
r = pm.Beta("r", alpha=alpha[0], beta=beta[0], shape=n, testval=r_obs)
mu_R = r * X_shared
R = pm.Poisson('R', mu=mu_R, observed=R_obs)
br = pm.Beta("br", alpha=alpha[1], beta=beta[1], shape=n, testval=br_obs)
# y is a factor of R and br
y_mu = R * br
y = pm.Poisson('y', mu=y_mu, observed=y_obs)
```

If I sample from this model’s posterior after replacing X_shared with new values, only variable `R`

, which is directly dependent on `X_shared`

will reflect the updated X values.

It seems like I *can* get this work if I replace `R`

in the definition of `y_mu`

with either `R_mu`

, which seems to make `y_mu`

directly dependent on `X_shared`

, with `r * X_shared`

explicitly.

However, in my actual model, not this toy model, this leads to much slower sampling rates. Presumably that’s b/c the sampling of `y`

is now dependent on add’l variables and can’t take advantage of `R`

.

Question: am I approaching this correctly? Any obvious errors in the above model? Or, is the model below the only way to to do? Any ideas much appreciated!!

```
with pm.Model() as model_predictable:
sd_alpha = pm.HalfNormal("sd_a", sd=1, shape=2)
sd_beta = pm.HalfNormal("sd_b", sd=1, shape=2)
alpha = pm.HalfNormal("alpha", sd=sd_alpha, shape=2)
beta = pm.HalfNormal("beta", sd=sd_beta, shape=2)
r = pm.Beta("r", alpha=alpha[0], beta=beta[0], shape=n, testval=r_obs)
R_mu = r * X_shared
R = pm.Poisson('R', mu=R_mu, observed=R_obs)
br = pm.Beta("br", alpha=alpha[1], beta=beta[1], shape=n, testval=br_obs)
y_mu = r * X_shared * br
y = pm.Poisson('y', mu=y_mu, observed=y_obs)
```