%posterior predictive - 
% P(X_n+1) 
p
post_mean = mean(posterior_sim);
frate(1:100) = post_mean(2) + post_mean(1)*(exp(-post_mean(3).*(1:100)));

%calculate the expected value of the next observation
post_prediction = expected_function(post_mean(1),post_mean(2),post_mean(3));
var_post_pred = variances_func(post_mean(1),post_mean(2),post_mean(3));

display('#failed in month I+1, #predicted, 2*SD')
[K(1,(I+1)) post_prediction(I+1) 2*sqrt(var_post_pred(I+1))]