#
# Genetic linkage example from Chapter 9, Exercises 3 and 8
#
niter <- 1
minitial <- 20
etadivs <- 1000
mag <- 10
scale <- 4
x <- c(461,130,161,515)
pagethrow <- 12
m <- minitial
eta <- runif(m)
y1 <- rbinom(m,x[1],eta/(eta+1))
y4 <- rbinom(m,x[4],(1-eta)/(2-eta))
for (t in 1:niter)
  {
     mold <- m
     if (t > 30)
       m <- 200
     if (t > 40)
       m <- 1000
     i0 <- floor(runif(m,0,mold))
     eta  <- rbeta(m,y1[i0]+x[2]+1,x[3]+y4[i0]+1)
     y1 <- rbinom(m,x[1],eta/(eta+1))
     y4 <- rbinom(m,x[4],(1-eta)/(2-eta))
  }
p <- rep(0,etadivs)
for (etastep in 1:(etadivs-1)){
   eta <- etastep/etadivs
   term <- exp((y1+x[2])*log(eta)+(x[3]+y4)*log(1-eta)
                 +lgamma(y1+x[2]+x[3]+y4+2)-lgamma(y1+x[2]+1)
                 -lgamma(x[3]+y4+1))
   p[etastep] <- p[etastep] + sum(term)/m
   if ((0.43 <= eta)&&(eta <= 0.45))
     cat(" p[",eta,"] = ",p[etastep],"\n");
}
cat("\n");
cat("Maximum at",order(p)[etadivs]/etadivs,"\n\n")
steps <- etadivs/mag
for (bigstep in 1:steps)
  {
     eta <- bigstep/steps
     if ((0.3 <= eta)&&(eta <= 0.55))
       {
          if (10*eta-floor(10*eta+0.001) < 0.001) cat(eta)
          else cat("   ")
          cat("|")
          for (k in 1:(scale*p[mag*bigstep])) 
            if (scale*p[mag*bigstep] >= 1) cat("*")
          cat("\n")
        }
  }