################################################################################################
## Code by Vanja Dukic, U of Colorado at Boulder
## Bayesian Statistics and Computing 2019
##
## Week 14 - Propagation of uncertainty 
## 
######################################################################3
# Say we have a sample of size N from the posterior
# from a logistic regression with a single predictor x
# Let's use the logistic regression from the homework problem 
# 
# Each sample point will give us one logistic curve
# Note for the logistic curve we only need alpha1 and alpha0

dim(par.thinned)
# 1112     4
# alpha0, alpha1, delta0, delta1


x=seq(-10,10,by=0.1)


alpha0.samples = par.thinned[,1]
alpha1.samples = par.thinned[,2]

theta.samples = matrix(0,length(alpha0.samples),length(x))

for (xi in 1:length(x))
{theta.samples[,xi] = exp(alpha0.samples + x[xi]*alpha1.samples)/(1+exp(alpha0.samples + x[xi]*alpha1.samples))
}

xmatrix = matrix(x,length(alpha0.samples),length(x),byrow=TRUE)

plot(x,theta.samples[1,],type='l', lty=3, cex=0.2,  col='grey', xlab='x',ylab='Probability')
par(new=TRUE)
plot(x,theta.samples[4,],type='l', lty=3, cex=0.2, col='grey',axis='')

matplot(t(xmatrix),t(theta.samples),type='l', lty=3, cex=0.2, col='grey',xlab='x',ylab='Probability',main='Posterior Samples')

###########################################################
# Compare this to the sample from the normal
require('MASS')
Opt.sample = mvrnorm(1000,OptTheta$par,solve(OptTheta$hessian))
Opt.alpha0.samples = Opt.sample[,1]
Opt.alpha1.samples = Opt.sample[,2]
Opt.theta.samples = matrix(0,length(Opt.alpha0.samples),length(x))
xmatrix = matrix(x,length(Opt.alpha0.samples),length(x),byrow=TRUE)

for (xi in 1:length(x))
{Opt.theta.samples[,xi] = exp(Opt.alpha0.samples + x[xi]*Opt.alpha1.samples)/(1+exp(Opt.alpha0.samples + x[xi]*Opt.alpha1.samples))
}

x11()
matplot(t(xmatrix),t(Opt.theta.samples),type='l', lty=3, cex=0.2, col='grey',xlab='x',ylab='Probability',main='Normal Approx')


x11()
par(mfrow=c(1,2))
hist(Opt.theta.samples[,which(x==3)],main='Normal Approx',xlab='Prob for x=3')
hist(theta.samples[,which(x==3)],main='Posterior Samples',xlab='Prob for x=3')

sum(Opt.theta.samples[,which(x==3)] > .5)/length(Opt.theta.samples[,which(x==3)])
sum(theta.samples[,which(x==3)] > .5)/length(theta.samples[,which(x==3)])