####################################
# Section 6.4 --turn in material described at the end 
#Numerical integration 
#  Example of generating a sample of size NN 
# from the joint posterior distribution of a normal model 
#  with uninformative prior. 
# Our data is going to be generated artificially, to 
# convince ourselves that this works 
####################################
yy=rnorm(1000,0,sqrt(4))   # this is our data 
  nn=length(yy)   # gives us sample size of the data
 ss=var(yy)       # gives us variance of the data
 ybar=mean(yy)    # gives us mean of the data
  AA=(nn-1)/2     # prepare the arguments of the inverse gamma 
 BB=0.5*(nn-1)*ss
########## Open space to put the generated mu and sigma squared ######
mu=rep(NA, 5000)           #   we will generate 5000 simulated mus
 sigmasqr=rep(NA,5000)   # 
####### We have two ways of generating the inverted gamma ############
library(MCMCpack)   # if you have installed MCMCpack you can use rinvgamma
 for(m in 1:5000){    # this loop will give the simulations from the joint posterior
 sigmasqr[m]=1/(rgamma(1,AA)/BB)    # This generates the same as next command
# sigmasqr[m]= rinvgamma(1, AA, BB)  # use this if you downloaded MCMCpack

#### For each sigma squared generate a mu  ################
 mu[m]=rnorm(1,ybar,sqrt(sigmasqr[m]/nn))   # plug in sigmasqr and simulate mu given it
 }
 
 ############## Plot marginals for mu and sigma squared ############

  par(mfrow=c(2,1))  # create space to put the plots. 
 hist(mu,density=-1,yaxt="n",yalb="",xlab="mu",cex=2,nclass=50)  
 hist(sqrt(sigmasqr),density=-1,yaxt="n",ylab="",xlab="sigma",cex=2,nclass=50)
 
 ############ Things to do on your own#####################
 #
 # (a) Copy paste the plots into your file.
#
# (b) Summarize the marginal posterior distribution for mu. Find marginal mean of mu, 
#      marginal variance, intervals,  some probabilities. 
#####################################################
 












######################################################
# 6.5.1 We generate the posterior distribution of the difference of means #
# ####################################################


muc= rt(1000,31)*(0.24/sqrt(32))+1.103   # this generates from the posterior for muc
mut= rt(1000,35)*(0.24/sqrt(36))+1.173   # this generates from the posterior for muc
mudif=muc-mut    # the difference between means

summary(mudif)     # will give us some summary statistics 
sd(mudif)   # gives the standard deviation 

layout(rbind(c(1,3), c(0,0), c(2,2)), heights=c(2,lcm(0.5),1),respect=TRUE)
hist(muc, xlab="muc", main="Posterior for mean of control group")
hist(mudif, xlab="muc-mut")
post.interval=quantile(mudif, c(0.025, 0.975))
abline(v=post.interval)
hist(mut, xlab="mut", main="Posterior for mean of treatment group")

prob.zero=(sum(mudif<0))/1000 
prob.zero




####### To do on your own #################
# Copy paste the histograms in your file. Add code to  summarize the 
# distribution of the posterior difference (mudif)
# between the means of the two groups (mean, variance). In addition to 
# Interpret the result obtained by prob.zero... what is the conclusion you 
# reach in this study? 
#######################################