You model will run into problem quite easily, as the model_err is not bounded: if model_err**2 > y_est*(1 - y_est) the Beta distribution is undefined. You can try doing:
model_err = pm.Uniform('model_err',lower=0, upper=y_est*(1-y_est), shape=(2481,))
eg, see linke below about infering a comment rate - with reparameterization using the mean http://stats.stackexchange.com/questions/67443/does-the-beta-distribution-have-a-conjugate-prior/67700
http://stats.stackexchange.com/questions/12232/calculating-the-parameters-of-a-beta-distribution-using-the-mean-and-variance/12239
You can also reparameterized it using the mode, see John K. Kruschke 2015, cf. Eqn 9.4, pp. 223, also pp. 236 Figure 9.7 for the estimation of kappa using a Gamma prior.
Otherwise, try the logit-normal model as @aseyboldt said in Beta Regression in pyMC3