# # 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) } }