########## R script: binomBetaChk ##########

# For checking MCMC analysis of a simple binomial/beta
# situation against the exact analytic answer.

library(BRugs)
library(KernSmooth)

allData <- list(n=20,y=c(0,0,1,0,1,0,0,0,0,0,0,1,0,0,0,0,0,0,0,0))

parInit <- list(list(p=mean(allData$y)))

BRugsFit(data=allData,inits=parInit,parametersToSave="p",nBurnin=1000,
         nIter=10000,modelFile="binomBetaModel.txt",numChains=1)

pSamp <- samplesSample("p")

# Plot the estimated posterior density estimated from the data.

estPostDest <- bkde(pSamp,bandwidth=dpik(pSamp))
par(mfrow=c(1,1))
plot(estPostDest,type="l",xlab="probab. of success",ylab="density",
     main="estimated posterior density",lwd=3,col="red",ylim=c(0,5.1))

xg <- estPostDest$x
exaPostDest <- dbeta(xg,4,18)
lines(xg,exaPostDest,col="blue")

legend(0.4,4,c("MCMC est.","exact"),lwd=c(3,1),col=c("red","blue"))

cat("\n\n Exact posterior mean is\n")
print(4/22)

cat("\n\n Sample mean of MCMC sample is\n")
print(mean(pSamp))


########## End of binomBetaChk ##########