rm(list = ls())
setwd("c:\\E-drive\\Books\\LOGLIN3\\BAYES\\")
load("post-c-samp.Rda")
# beta.c

# Range of predictions
prediction_temps <-  (300:800)/10

# Center prediction times
tau=c(53,57,58,63,66,67,67,67,68,69,70,70,70,70,72,
73,75,75,76,76,78,79,81)
temp_center=prediction_temps - mean(tau)
# Prediction Model Matrix
# Create a vector of 1s
J=1+temp_center - temp_center
# Create matrix X_f
Xf=c(J,temp_center)
Xf=matrix(Xf,ncol=2)
Xf

# compute posterior mean predictions
Pr=Xf %*% beta.c
Pr=exp(Pr)/(1+exp(Pr))
MP=rowMeans(Pr)
# look at prediction probabilities
matrix(c(prediction_temps,MP),ncol=2)

# posterior predictive probabilities:
#         5%, 95% quantiles and median
prediction_quantiles <- apply(Pr, 1,
                     quantile, probs=c(0.05,0.5,0.95))

# Now prepared to do Fig 13.8 because it is
# all based on the basic sample of 10000 betas
# Easy to replace posterior means
# (MP) with posterior medians (Median)
# Fig 13.8, Fletch
predictions <- data.frame( prediction_temps, t(prediction_quantiles), MP )
names(predictions) <- c("Temp", "Lower", "Median", "Upper","MP")
library(ggplot2)
ggplot( predictions ) +
  geom_line( aes(x=Temp, y=Lower),
  color="#aa0000", linetype=2 ) +
  geom_line( aes(x=Temp, y=MP),
  color="#aa0000", linetype=1 ) +
  geom_line( aes(x=Temp, y=Upper),
  color="#aa0000", linetype=2 ) +
  coord_cartesian( ylim=c(0,1) ) +
  theme_bw() +
  xlab("Temperature") +
  ylab("Failure Probability") +
  ggtitle( expression( paste(
  "O-Ring failure probability by temperature" ) ) )

# Additional computations for Fig 13.7

# Case 18 deleted, mean predictions
load("post-c-samp-18.Rda")
Pr=Xf %*% beta.c
Pr=exp(Pr)/(1+exp(Pr))
MP_18=rowMeans(Pr)

# MLE predictions from Chap 2
beta1mle=15.0429
beta2mle=-.2322
#uncentered MLEs
lp= beta1mle + prediction_temps*beta2mle
mle_predictions <- exp(lp) / ( 1 + exp(lp) )

# Fig 13.7
plot(prediction_temps,mle_predictions,type="l",
     xlim=c(30,85),
     ylab="Probability",xlab="Temperature",lty=3)
lines(prediction_temps,MP,type="l",lty=1)
lines(prediction_temps,MP_18,type="l",lty=2)
legend("bottomleft",
     c("MLE","Bayes: full data","Bayes: no case 18"),
     lty=c(3,1,2))