#
# Chained data augmentation - Example from Casella and George
#
nr <- 50
m <- 500
k <- 10
n <- 16
alpha <- 2.0
beta <- 4.0
lambda <- 16.0
maxn <- 24
betabinomial <- function(x,n,alpha,beta)
{
   y <- log(choose(n,x))
   y <- y + lgamma(alpha + beta) - lgamma(alpha) - lgamma(beta)
   y <- y + lgamma(x + alpha) + lgamma(n - x + beta) - 
           lgamma(alpha + beta + n)
   y <- exp(y)
   return(y)
}
cat("\n")
cat("Based on 'Explaining the Gibbs sampler', C. Casella      \n")
cat("and E.I. George, Amer. Statist. 46 (3) (1992), 167-174. \n")
h <- rep(0,n+1)
fe <- rep(0,n+1)
for (i in 1:m)
      {
         y <- runif(1);
         for (j in 1:k)
            {
               x <- rbinom(1,n,y)
               newalpha <- x + alpha
               newbeta <- n - x + beta
               y <- rbeta(1,newalpha,newbeta)
            }
         for (t in 0:n)
            {
               if (t == x)
                  h[t+1] <- h[t+1] + 1
               term <- choose(n,t)*exp(t*log(y)+(n-t)*log(1-y))
               fe[t+1] <- fe[t+1] + term
            }
      }
cat("\n")
cat("Histogram (cf. Fig. 1)) \n")
cat("         t       Obs      Exp     Diff       Ratio         Comp of X2 \n")
cat("\n")
x2h <- 0
bbe <- rep(0,n+1)
bb <- rep(0,n+1)
for (t in 0:n)
      {
         bbe[t+1] <- m*betabinomial(t,n,alpha,beta)
         bb[t+1] <- round(bbe[t+1])
         diff <- h[t+1] - bb[t+1]
         ratio <- h[t+1]/bbe[t+1]
         compx2 <- (h[t+1]-bbe[t+1])*(h[t+1]-bbe[t+1])/bbe[t+1]
         x2h <- x2h + compx2
         if (t < 10) cat(" ")
         cat("        ",t,"     ")
         if (h[t+1] < 10) cat(" ")
         cat(h[t+1],"      ")
         if (bb[t+1]<10) cat(" ")
         cat(bb[t+1],"     ")
         if (diff >= 0) cat(" ")
         if (abs(diff) < 10) cat(" ")
         cat(diff,"      ",ratio,"    ",compx2,"\n")
      }
cat("\n")
cat("Chi-squared equals",x2h,"on",n,"degrees of freedom \n")
cat("\n")
cat("Estimated densities (cf. Fig. 3) \n")
cat("\n")
cat("         t       Obs      Exp     Diff       Ratio         Comp of X2 \n")
cat("\n")
x2f <- 0
f <- rep(0,n)
for (t in 1:n)
      {
         f[t+1] <- round(fe[t+1])
         diff <- f[t+1] - bb[t+1]
         ratio <- f[t+1]/bbe[t+1]
         compx2 <- (f[t+1]-bbe[t+1])*(f[t+1]-bbe[t+1])/bbe[t+1]
         x2f <- x2f + compx2
         if (t < 10) cat(" ")
         cat("        ",t,"     ")
         if (f[t+1] < 10) cat(" ")
         cat(f[t+1],"      ")
         if (bb[t+1]<10) cat(" ")
         cat(bb[t+1],"     ")
         if (diff >= 0) cat(" ")
         if (abs(diff) < 10) cat(" ")
         cat(diff,"      ",ratio,"    ",compx2,"\n")
      }
cat("\n")
cat("Chi-squared equals",x2f,"on",n,"degrees of freedom. \n")
hp <- rep(0,(maxn+1))
fep <- rep(0,(maxn+1))
for (i in 1:m)
      {
         y <- 0.5
         nn <- (1-y)*lambda;
         for (j in 1:k)
            {
               x <- rbinom(1,nn,y)
               newalpha <- x + alpha
               newbeta <- nn - x + beta
               y <- rbeta(1,newalpha,newbeta)
               nn <- x + rpois(1,(1-y)*lambda)
            }
         for (t in 0:maxn)
            {
               if (t == x)
                  hp[t+1] <- hp[t+1] + 1
               if (t <= nn)
                  {
                      term <- choose(nn,t)*exp(t*log(y)+
                             (nn-t)*log(1-y))
                      fep[t+1] <- fep[t+1] + term
                  }
            }
      }
cat("\n\n")
cat("Histogram (n random) \n")
cat("\n")
cat(" t       Obs    Histogram \n")
cat("\n")
practmaxn <- 4*n/3
for (t in 0:(practmaxn+1))
      {
         if (t < 10) cat(" ")
         cat(t,"      ")
         if (hp[t+1] < 10) cat(" ")
         cat(hp[t+1],"    ")
         if (hp[t+1] > 0)
           for (j in 1:hp[t+1]) cat("*")
         cat("\n")
      }
cat("\n")
pp <- 0:(length(hp)-1)
barplot(hp,ylim=c(0,80),names=pp,main="Histogram (n random)")
muh <- sum(hp*pp)/sum(hp)
ex2h <- sum(hp*pp*pp)/sum(hp)
sigma2h <- ex2h - muh*muh
sigmah <- sqrt(sigma2h)
cat("\nMean",muh,"s.d.",sigmah,"\n\n")
cat("Estimated densities (n random; cf. Fig. 5) \n")
cat("\n")
cat(" t       Obs    Estimate \n")
cat("\n")
x2f <- 0
fp <- rep(0,practmaxn)
for (t in 1:practmaxn)
      {
         fp[t+1] <- round(fep[t+1])
         if (t < 10) cat(" ")
         cat(t)
         cat("       ")
         if (fp[t+1] < 10) cat(" ")
         cat(fp[t+1],"    ")
         if (fp[t+1] > 0)
           for (j in 1:fp[t+1]) cat("*")
         cat("\n")
      }
p <- 0:(length(fp)-1)
barplot(fp,ylim=c(0,60),names=p,
        main="Estimated densities (n random; cf. Fig. 5)")
mu <- sum(fp*p)/sum(fp)
ex2 <- sum(fp*p*p)/sum(fp)
sigma2 <- ex2 - mu*mu
sigma <- sqrt(sigma2)
cat("\nMean",mu,"s.d.",sigma,"\n\n")