rm(list = ls())
#Read posterior samples of beta
# Assumes they are in a p x t matrix
setwd("c:\\E-drive\\Books\\LOGLIN3\\BAYES\\")
load("trauma-samp-thin.Rda")
#load("trauma_samp.Rda")
beta
# establish t (number of iterates)
iterates=length(beta[1,])

# define prior hyperparameters
a=c(1.1,3,5.9,1.3,1.1,1.5)
b=c(8.5,11,1.7,12,4.9,5.5)

#Enter Xtilde
Xtp=c(
 1, 25, 7.84, 60, 0,  0,
 1, 25, 3.34, 10, 0,  0,
 1, 41, 3.34, 60, 1, 60,
 1, 41, 7.84, 10, 1, 10,
 1, 33, 5.74, 35, 0,  0,
 1, 33, 5.74, 35, 1, 35)
Xtilde=t(matrix(Xtp,6,6))
Xtilde

# generate posterior samples of ptildes
# Compute posterior linear predictors and probabilities
ptilde_posterior= Xtilde %*% beta
ptilde_posterior=exp(ptilde_posterior)/(1+exp(ptilde_posterior))

# Matrix for plots
par(mfrow=c(2,3))

# plot
# generate prior sample within 2nd call of density
 for(i in 1:6){
post_p=density(ptilde_posterior[i,])
prior_p=density(rbeta(iterates,a[i],b[i]) )
plot(prior_p,main=c(i),xlab="ptilde", ylim=c(0,10), xlim=c(0,1) )
lines(post_p,lty=5)
legend("topright",c("prior","post"),lty=c(1,5))
}