library(MCMCpack) library(MASS) hist(newcomb,breaks=30) y <- newcomb n <- length(y) #Estimating the model and drawing n.sims posterior simulations #We assume the standard uninformative prior; see Section 9.2 in lecture notes n.sims <- 1000 sigma2.sim <- (n-1)*var(y)/rchisq(n.sims,df=n-1) mu.sim <- rnorm(n.sims,mean(y),sigma2.sim/n) #Generating n.sims replicates of the original data set #and computing a test quantity and a test statistics from each n <- length(y) y.rep <- matrix(NA,n.sims,n) test.quant <- numeric(n.sims) test.quant.rep <- numeric(n.sims) test.stat.rep <- numeric(n.sims) test.stat <- min(y) for(i in 1:n.sims){ y.rep[i,]<- rnorm(n,mu.sim[i],sqrt(sigma2.sim[i])) test.quant.rep[i] <- mean((y.rep[i,]-mu.sim[i])^2) test.quant[i] <- mean((y-mu.sim[i])^2) test.stat.rep[i] <- min(y.rep[i,]) } #posterior predicive p-values: mean(test.quant.rep>=test.quant) #[1] 0.506 mean(test.stat.rep>=test.stat) #[1] 1 #Graphs #Replicated data sets op <- par(mfrow=c(5,4)) par(plt=c(0.2,0.8,0.2,0.8)) for(i in 1:20){ hist(y.rep[i,],main="",ylab="") } par(op) #Replicated test statistic (min) and test quantity (var) op <- par(mfrow=c(1,2)) hist(test.stat.rep,xlim=c(-50,20),main="") lines(c(test.stat,test.stat),c(0,350),lwd=2,col=2) plot(test.quant,test.quant.rep, xlim=c(0,250),ylim=c(0,250)) abline(a=0,b=1) par(op)