example1 <- function(x=10, n=28, a=1, b=1, N=5000){

theta=seq(0.001, 0.999, l=200)
theta.prior=dbeta(theta, a, b)
theta.posterior=dbeta(theta, a+x, b+n-x)

par(mfrow=c(2,1),cex=0.8)
plot(theta, theta.posterior, ylab="Density", type="l", col="blue")
lines(theta, theta.prior, type="l", col="red", lty=2)
legend(0.8,3,legend=c("posterior","prior"),lty=1:2)

theta=rbeta(N, a+x, b+n-x)
eta=log((1-theta)/theta)

hist(eta,probability=T,col="lightblue",border="darkblue")

cat("Posterior summary of eta:",fill=T)
print(summary(eta))
cat(c("posteriot sd of eta:",round(sd(eta),4)),fill=T)
cat(c("MC se of eta:",round(sd(eta)/sqrt(N),4)),fill=T)

cat("Posterior percentiles of eta:",fill=T)
print(round(quantile(eta,c(1:9)/10),4))

}