#
# Genetic linkage example from Sections 9.2 and 9.3
#
niter <- 50
minitial <- 20
etadivs <- 1000
mag <- 10
scale <- 8
x <- c(125,18,20,34)
pagethrow <- 12
m <- minitial
eta <- runif(m)
y <- rbinom(m,x[1],eta/(eta+2))
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,y[i0]+x[4]+1,x[2]+x[3]+1)
     y <- rbinom(m,x[1],eta/(eta+2))
  }
p <- rep(0,etadivs)
for (etastep in 1:(etadivs-1)){
   eta <- etastep/etadivs
   term <- exp((y+x[4])*log(eta) + (x[2]+x[3])*log(1-eta)
                 +lgamma(y+x[2]+x[3]+x[4]+2)-lgamma(y+x[4]+1)
                 -lgamma(x[2]+x[3]+1))
   p[etastep] <- p[etastep] + sum(term)/m
   if ((0.62 <= eta)&&(eta <= 0.635))
     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.45 <= eta)&&(eta <= 0.81))
       {
          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")
        }
}
vals <- 1:etadivs/etadivs
estmean <- sum(p*vals)/etadivs
estss <- sum(p*vals*vals)/etadivs
estsd <- sqrt(estss - estmean*estmean)
cat("Estimated mean",estmean,"\n")
cat("Estimated s.d.",estsd,"\n")