#'simulate examples for WAPL
#'This is the generative function for WAPL
#'@param scenario: takes value 1-4, represent different scenarios
#'@param N: number of subjects
#'@param sigma: value of noise
#'@export
#'@return a list contains generated data and related information
gSim<-function(scenario=c(1,2,3,4), N=200, sigma=0.2){
#p=50 predictors in 20 blocks of 3; the later 40 have all coef of zero
#within group the data are uniform with correlation 0.2
A=sample(c(1,-1), N, replace=T, prob=c(0.5,0.5))
if(scenario==1){
X = Reduce(cbind, lapply(1:50, function(...){
##a block of uniform variable, within group correlation is 0.2
matrix(runif(N*2, 0, 1), nrow=N) + matrix(rep(runif(N, 0, 1), 2), nrow=N, byrow=F)
}))
colnames(X)=paste("X",c(1:100),sep="")
Ym= 2*X[,3]+ 2*pi*(sin(pi*X[,4]))
Yc= 3*pi*(sin(pi*X[,1]*2)) - 0.5*(X[,2]-0.5)^2
}else if(scenario==2){
X = Reduce(cbind, lapply(1:10, function(...){
##a block of uniform variable, within group correlation is 0.2
matrix(runif(N*2, 0, 1), nrow=N) + matrix(rep(runif(N, 0, 1), 2), nrow=N, byrow=F)
}))
colnames(X)=paste("X",c(1:20),sep="")
Ym= 2*X[,3]+ 2*pi*(sin(pi*X[,4]))
Yc= 3*pi*(sin(pi*X[,1]*2)) - 0.5*(X[,2]-0.5)^2
}else if(scenario==3){
X = Reduce(cbind, lapply(1:50, function(...){
##a block of uniform variable, within group correlation is 0.2
matrix(runif(N*2, 0, 1), nrow=N) + matrix(rep(runif(N, 0, 1), 2), nrow=N, byrow=F)
}))
colnames(X)=paste("X",c(1:100),sep="")
Ym= 4*X[,5]+ 4*pi*(sin(pi*X[,6]))
Yc= 5*pi*(sin(pi*X[,1]*2)) - 1*(X[,2]-0.5)^2 + 10*(sin(pi*X[,3])) + 2*cos(pi*X[,4]*2)
} else if(scenario==4){
X = Reduce(cbind, lapply(1:10, function(...){
##a block of uniform variable, within group correlation is 0.2
matrix(runif(N*2, 0, 1), nrow=N) + matrix(rep(runif(N, 0, 1), 2), nrow=N, byrow=F)
}))
colnames(X)=paste("X",c(1:20),sep="")
Ym= 4*X[,5]+ 4*pi*(sin(pi*X[,6]))
Yc= 5*pi*(sin(pi*X[,1]*2)) - 1*(X[,2]-0.5)^2 + 10*(sin(pi*X[,3])) + 2*cos(pi*X[,4]*2)
} else {break;}
noise=rnorm(N, sd=sigma)
s2n=var(Yc) /(var(Ym)) #+ var(noise))
Y=Ym + A*Yc + noise
optA={Yc > 0}*2-1
Y0=Ym - Yc #Everyone follows A=-1
Y1=Ym + Yc #Everyone follows A=1
optY=Ym + abs(Yc)
tmp=table(optA)
out=list("X"=X, "A"=A, "Y"=Y, "optA"=optA, "Y1"=Y1, "Y0"=Y0, "optY"=optY,
"s2n"=s2n, "PrOptA1"= tmp[2]/N, "EY"=mean(optY), "EYopt0"=mean(optY[optA==-1]),
"EYopt1"=mean(optY[optA==1]), "EYAll0"=mean(Y0), "EYAll1"=mean(Y1),
"MedY"=median(optY), "MedYopt0"=median(optY[optA==-1]),
"MedYopt1"=median(optY[optA==1]), "MedYAll0"=median(Y0), "MedYAll1"=median(Y1))
#return(unlist(out[8:17]))
return(out)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.