Thanks for the example. But what if the X and y are already modeled?
X = pm.Normal('X', mu = hyper_x_mu, 1, shape=(10))
y = pm.Normal('y', mu = hyper_y_mu, sigma = hyper_y_sigma, shape=(10))
Where hyper_x_mu, hyper_y_mu and hyper_y_sigma are RVs from upstream models.
Should I do this… or?
residuals = pm.Normal('residuals', mu=0, sigma=100, shape=(10))
alpha = pm.Normal('alpha', 0, 1)
beta = pm.Normal('beta', 0, 1)
should_be_zero = ((y + residuals) - alpha) / beta - 1
#use potential to constrain 'should_be_zero' to zero?