Computing rmse with y_rep matrix

I have made a simple bayesian negative binomial model for modelling count data based on this tutorial: http://www.stat.columbia.edu/~gelman/bda.course/_book/a-case-study-in-bayesian-workflow.html. Now, I have my generated in sample y_rep and real y and when I make a ppc_dens_overlay plot the generated model looks like to fit my data pretty well (see attachment). However, when I calculate the rmse between the generated y_rep and y, I get a very high value. So, I did something wrong. To calculate the rmse, I first take the column means of my y_rep matrix, but since I have lots of zeros but also large values the estimated mean is far of the real y value. Does anyone know how to correctly compute the rmse using the y_rep matrix and to compute point forecasts in this case?

plot_mc.pdf (447.6 KB)