# This function is one atomic simulation
# of the normal setting
# which returns the following estimates for AUC:
# thin version does not include:
# SIMEX "naive"
# SIMEX on averages
# SIMEX on EBLUPs
atomic.sim.thin <- function(m,n,p,muX,muY,
sigmaX,sigmaY,sigmaEps,sigmaEta,
K,LAM1,LAM2)
{
X <- rnorm(m,muX,sigmaX)
Y <- rnorm(n,muY,sigmaY)
# additive normal ME's with var sigmaEps, sigmaEta
eps <- matrix(rnorm(n = m*p,mean = 0, sd = sigmaEps),nrow = m,ncol = p)
eta <- matrix(rnorm(n = n*p,mean = 0, sd = sigmaEta),nrow = n,ncol = p)
x <- X + eps
y <- Y + eta
# eblup
eblup.est <- eblup.auc.REML(x,y)
e <- eblup.est[[1]] # the MW estimate from EBLUPS
sEps <- eblup.est[[2]][1] # REML estimate of sigma_eps
sEta <- eblup.est[[2]][2] # REML estimate of sigma_eta
ebX <- eblup.est$blups1
ebY <- eblup.est$blups2
# simex
#simex <- simex.auc.naive.RM(y,x,K=K,LAM1 = LAM1,LAM2 = LAM2,tau1.reml = sEta,tau2.reml = sEps)
#s1 <- simex[1] # linear
#s2 <- simex[2] # quadratic
# Reiser
#r <- reiser.2000.REML(x,y,sX,sY)
# Simex on averages
simex.aver <- simex.auc.reml(rowMeans(x),rowMeans(y),K = K,LAM1 = LAM1,LAM2 = LAM2,sigmaX = sEps,sigmaY = sEta)
sa1 <- simex.aver[1]
sa2 <- simex.aver[2]
# simex on eblups
simex.eblups <- simex.auc.reml(ebX,ebY,K = K,LAM1 = LAM1,LAM2 = LAM2,sigmaX = sEps,sigmaY = sEta)
vec <- c(e,sa1,sa2,simex.eblups)
names(vec) <- c("Eblup","average1","average2","simex-eblup1","simex-eblup2")
return(vec)
}
# Experiment:
# m <- n <- 200
# p <- 3
# muX <- 1.8
# muY <- 0
# sigmaX <- sigmaY <- 1
# sigmaEps <- sigmaEta <- 1
# K <- 1000
# LAM1 <- LAM2 <- seq(0.1,1.9,by=0.2)
#
# a <- Sys.time()
# atomic.sim.thin(m=m,n=n,p=p,muX=muX,muY=muY,sigmaX = sigmaX,sigmaY = sigmaY,sigmaEps = sigmaEps,sigmaEta = sigmaEta,
# K=10,LAM1 = LAM1,LAM2=LAM2)
# Sys.time()-a
#
# microbenchmark::microbenchmark(atomic.sim.thin(m=m,n=n,p=p,muX=muX,muY=muY,sigmaX = sigmaX,sigmaY = sigmaY,sigmaEps = sigmaEps,sigmaEta = sigmaEta,
# K=1000,LAM1 = LAM1,LAM2=LAM2),times = 5)
#
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.