# NOT RUN {
d <- 4
n <- 200
sigma <- 0.5
delta <- 1
height <-1
arange <- c(0,5)
triangle <- function(a,height){
y <- exp(-a^2/((1/2)^2))*height
return(y)
}
mu.mod<-function(a,l,delta,height){
mu <- as.numeric(l%*%c(0.2,0.2,0.3,-0.1*delta))+
triangle(a-2.5,height)+a*(-0.1*l[,1]+0.1*delta*l[,4])
return(mu)
}
l <- matrix(rnorm(n*d),ncol=d)
l[,4] <- ifelse(l[,4]>0,1,0)
colnames(l) <- paste("l",1:4,sep="")
logit.lambda <- as.numeric(l%*%c(0.1,0.1,-0.1,0))
lambda <- exp(logit.lambda)/(1+exp(logit.lambda))
a <- rbeta(n, shape1 = lambda, shape2 =1-lambda)*5
mu <- mu.mod(a,l,delta,height)
residual.list <- rnorm(n,mean=0,sd =sigma)
y <- mu+residual.list
class_label <- l[,4]
ylist <- split(y,class_label)
alist <- split(a,class_label)
pilist <- split(pmin(dbeta(a/5,shape1=lambda,shape2=1-lambda)/5,100),class_label)
mulist <- split(mu,class_label)
varpilist <- list()
malist <- list()
for(c in c(0,1)){
ac <- a[class_label==c]
lc <- l[class_label==c,]
logit.lambdac <- as.numeric(lc[rep(1:nrow(lc),nrow(lc)),]%*%c(0.1,0.1,-0.1,0))
lambdac <- exp(logit.lambdac)/(1+exp(logit.lambdac))
varpic <- colMeans(matrix(pmin(dbeta(rep(ac,each=length(ac))/5,
shape1=lambdac,
shape2 = 1-lambdac)/5,100),nrow=length(ac)))
mac <- colMeans(matrix(mu.mod(rep(ac,each=length(ac)),
lc[rep(1:nrow(lc),nrow(lc)),],
delta,height),
nrow=length(ac)))
varpilist[[as.character(c)]]<-varpic
malist[[as.character(c)]] <- mac
}
out <- drdrtest_em.base(ylist,alist,pilist,varpilist,mulist,malist,arange)
# }
Run the code above in your browser using DataLab