%posterior predictive - 
% P(X_n+1) 
p
for i=1:MM
    post_pred_sample(i,:) = expected_function(posterior_sim(i,1),posterior_sim(i,2),posterior_sim(i,3));
    post_var_sample(i,:) = variances_func(posterior_sim(i,1),posterior_sim(i,2),posterior_sim(i,3));
end
post_mean = mean(post_pred_sample);
post_v = mean(post_var_sample);

post_prediction(count) = post_mean(I+1);
post_var(count) = post_v(I+1)
upper(count) = post_prediction(count) + 2*sqrt(post_var(count));
lower(count) = post_prediction(count) - 2*sqrt(post_var(count));
display('#failed in month I+1, #predicted, 2*SD')
[obs(1,(I+1)) post_prediction(count) 2*sqrt(post_var)]



