# We use the following functions:

OU.sim <- function(beta,markers,n.iter=10^4)
{
    require(MASS)
    V <- outer(markers,markers,function(x,y) exp(-beta*abs(x-y)))
    mu <- rep(0,length(markers))
    Z <- mvrnorm(n.iter,mu,V)
    return(Z)
}

add.qtl <- function(Z,beta,markers,q,xi)
{
    d <- dim(Z)
    if (length(markers) != d[2])
        stop("Number of columns of simulated matrix
        does not match the number of markers")
    mu <- xi*exp(-beta*abs(markers-q))
    Z <- sweep(Z,2,mu,"+")
    return(Z)
}


# Determine the threshold

markers <- seq(0,80,by=10)
beta.a <- 0.02
beta.d <- 0.04
Z.a <- Z.d <- NULL
for (i in 1:20) 
{
   Z.a <- cbind(Z.a,OU.sim(beta.a,markers))
   Z.d <- cbind(Z.d,OU.sim(beta.d,markers))
}
U.max <- apply(Z.a^2 + Z.d^2,1,max)
thresh <- quantile(U.max,0.95)


# Find the power
n <- 100
al <- 2
del <- 1
sig <- 1
q <- 35
xi.a <- sqrt(n/2)*al/sig
xi.d <- sqrt(n/4)*del/sig
Z.a <- OU.sim(beta.a,markers)
Z.d <- OU.sim(beta.d,markers)
Z.a1 <- add.qtl(Z.a,beta.a,markers,q,xi.a)
Z.d1 <- add.qtl(Z.d,beta.d,markers,q,xi.d)
mean(apply(Z.a1^2+Z.d1^2,1,max)> thresh)