#
# linkage.r (Genetic Linkage Example)
#

# Start algorithm
y <- c(125, 18, 20, 34)
# input data y
integral <- .2357695165e29
# The integral was calculated beforehand using MAPLE
# integral:= int((2+x)^125*(1-x)^(18+20)*x^34, x=0..1);
true <- function(x)((2+x)^y[1]*(1-x)^(y[2]+y[3])*x^y[4]/integral)
# calculate the true posterior
m <- 1600
# take 1600 samples
theta <- runif(m, min=0, max=1)
# get samples from prior distribution
x2 <- rbinom(m, y[1], theta/(theta+2))
# augmentation
posterior <- function(x)(1/m*sum(dbeta(x, x2+y[4], y[2]+y[3])))
# estimate posterior distribution
datax <- (0:1000)
for (i in 1:1001) datax[i] <- (i-1)/1000
datay <- (0:1000)
for (i in 1:1001) datay[i] <- posterior((i-1)/1000)
x11(record=T)
plot(datax, datay, ylim=c(0,8), type="l",
     xlab="theta", ylab="density", main=1)
# plot the estimate
curve(true, 0, 1, n=1001, add=TRUE)
# plot the true posterior to compare it to the estimate

for (j in 1:49)
# execute the algorithm 49 more times
{
    decomp <- sample(m, m, replace=TRUE)
    # decomposite the estimate according to [KG] section 6.4.2
    theta <- (0:1000)
    for (i in 1:1001) theta[i] <-
         rbeta(1, x2[decomp[i]]+y[4], y[2]+y[3])
    x2 <- rbinom(m, y[1], theta/(theta+2))
    if ((j==9)|(j==24)|(j==49))
    {
        posterior <-
             function(x)(1/m*sum(dbeta(x, x2+y[4], y[2]+y[3])))
        for (i in 1:1001) datay[i] <- posterior((i-1)/1000)
        x11(record=T)
        plot(datax, datay, ylim=c(0,8), type="l",
             xlab="theta", ylab="density", main=j+1)
        curve(true, 0, 1, n=1001, add=TRUE)
    }
}