## auxiliary functions for dense setting
# (01) scenario.2014CLX : two-means test (4 means, 8 covariances)
# (02) scenario.2013CLX : two-covariance test (4 covariances)
# (03) gen.cov.2014CLX : 8 covariance matrix generator (scenario.2014CLX)
# (04) gen.2013CLX.mine : sparse covariance difference setting
# (05) gen.2013CLX.dense : dense covariance difference setting
# (01) scenario.2014CLX ---------------------------------------------------
#' @keywords internal
#' @noRd
scenario.2014CLX <- function(n, p, no.mean=1, no.cov=1, htype=c("H0","H1")){
# put all the relevant functions here for clarity
scenario.2014CLX.mean1 = function(n, p){
m = max(1,floor(0.05*p))
vals = sample(c(-1,1), m, replace=TRUE)*sqrt(log(p)/n)
mu0 = rep(0,p)
mu0[sample(1:p, m, replace=FALSE)] = vals
return(mu0)
}
scenario.2014CLX.mean2 = function(n, p){
m = max(1,floor(0.05*p))
bound = sqrt(8*log(p)/n)
vals = runif(m, min=-bound, max=bound)
mu0 = rep(0,p)
mu0[sample(1:p, m, replace=FALSE)] = vals
return(mu0)
}
scenario.2014CLX.mean3 = function(n, p){
m = max(1,sqrt(p))
vals = sample(c(-1,1), m, replace=TRUE)*sqrt(log(p)/n)
mu0 = rep(0,p)
mu0[sample(1:p, m, replace=FALSE)] = vals
return(mu0)
}
scenario.2014CLX.mean4 = function(n, p){
m = max(1,sqrt(p))
bound = sqrt(8*log(p)/n)
vals = runif(m, min=-bound, max=bound)
mu0 = rep(0,p)
mu0[sample(1:p, m, replace=FALSE)] = vals
return(mu0)
}
scenario.2014CLX.cov1 = function(n, p){
Sig = array(0,c(p,p))
for (k in 1:floor(p/2)){
lowk = as.integer(2*(k-1)+1)
uppk = as.integer(2*k)
for (i in lowk:uppk){
for (j in lowk:uppk){
Sig[i,j] = 0.8
}
}
}
diag(Sig) = 1
return(Sig)
}
scenario.2014CLX.cov2 = function(n, p){
Sig = array(0,c(p,p))
for (i in 1:p){
for (j in 1:p){
Sig[i,j] = 0.6^(abs(i-j))
}
}
return(Sig)
}
scenario.2014CLX.cov3 = function(n, p){
Omega = array(0,c(p,p))
for (i in 1:(p-1)){
Omega[i,i+1] = 0.8
}
for (i in 1:(p-2)){
Omega[i,i+2] = 0.4
}
for (i in 1:(p-3)){
Omega[i,i+3] = 0.2
}
Omega = Omega + t(Omega)
diag(Omega) = 2
Sig = solve(Omega, diag(p))
return(Sig)
}
scenario.2014CLX.cov4 = function(n, p){
Omega = array(0,c(p,p))
for (i in 1:(p-1)){
for (j in (i+1):p){
theval = 0.6^(abs(i-j))
Omega[i,j] = theval
Omega[j,i] = theval
}
}
diag(Omega) = 1
Dhalf = diag(sqrt(runif(p, min=1, max=3)))
Sig = Dhalf%*%solve(Omega,Dhalf)
return(Sig)
}
scenario.2014CLX.cov5 = function(n, p){
Dhalf = diag(sqrt(runif(p, min=1, max=3)))
Ohalf = array(0,c(p,p))
for (k in 1:floor(p/2)){
lowk = as.integer(2*(k-1)+1)
uppk = as.integer(2*k)
for (i in lowk:uppk){
for (j in lowk:uppk){
Ohalf[i,j] = 0.8
}
}
}
diag(Ohalf) = 1
Omega = Dhalf%*%Ohalf%*%Ohalf%*%Dhalf
Sig = solve(Omega, diag(p))
}
scenario.2014CLX.cov6 = function(n, p){
SigStar = array(0,c(p,p))
for (k in 1:floor(p/2)){
lowk = as.integer(2*(k-1)+1)
uppk = as.integer(2*k)
for (i in lowk:uppk){
for (j in lowk:uppk){
SigStar[i,j] = 0.8
}
}
}
diag(SigStar) = 1
Dhalf = diag(sqrt(runif(p, min=1, max=3)))
E = array(0,c(p,p))
noff = sum(1:(p-1))
E[lower.tri(E)] = rbinom(noff, 1, 0.3)*runif(noff, min=-0.2, max=0.2)
E = E + t(E)
triplet = Dhalf%*%SigStar%*%Dhalf
delta = abs(min(eigen(triplet + E)$values))
Sig = triplet + E + delta*diag(p)
return(Sig)
}
scenario.2014CLX.cov7 = function(n, p){
SigStar = array(0,c(p,p))
for (i in 1:(p-1)){
for (j in (i+1):p){
theval = ((abs(i-j))^{-5})/2
SigStar[i,j] = theval
SigStar[j,i] = theval
}
}
diag(SigStar) = 1
Dhalf = diag(sqrt(runif(p, min=1, max=3)))
Sig = Dhalf%*%SigStar%*%Dhalf
return(Sig)
}
scenario.2014CLX.cov8 = function(n, p){
Dhalf = diag(sqrt(runif(p, min=1, max=3)))
Fmat = diag(p)
for (i in 1:(p-1)){
j = i+1
Fmat[i,j] = 0.5
Fmat[j,i] = 0.5
}
UmatOriginal = matrix(rnorm(p*p), nrow=p)
UmatOriginal = UmatOriginal%*%t(UmatOriginal)
U3 = qr.Q(qr(UmatOriginal))[,1:3]
Sig = Dhalf%*%(Fmat + (U3%*%t(U3)))%*%Dhalf
return(Sig)
}
# 1. generate Ha's mean
no.mean = as.integer(no.mean)
if (!(no.mean %in% c(1,2,3,4))){
stop("* scenario.2014CLX : 'no.mean' is from 1 to 4.")
}
mu.alpha = switch(no.mean,
"1" = scenario.2014CLX.mean1(n,p),
"2" = scenario.2014CLX.mean2(n,p),
"3" = scenario.2014CLX.mean3(n,p),
"4" = scenario.2014CLX.mean4(n,p))
# 2. generate common covariance
no.cov = as.integer(no.cov)
if (!(no.cov %in% 1:8)){
stop("* scenario.2014CLX : 'no.cov' is from 1 to 8.")
}
sig.alpha = switch(no.cov,
"1" = scenario.2014CLX.cov1(n,p),
"2" = scenario.2014CLX.cov2(n,p),
"3" = scenario.2014CLX.cov3(n,p),
"4" = scenario.2014CLX.cov4(n,p),
"5" = scenario.2014CLX.cov5(n,p),
"6" = scenario.2014CLX.cov6(n,p),
"7" = scenario.2014CLX.cov7(n,p),
"8" = scenario.2014CLX.cov8(n,p))
# generation and return
htype=match.arg(htype)
output = list()
if (htype=="H0"){
output$X = mvtnorm::rmvnorm(n, mean=rep(0,p), sigma=sig.alpha)
output$Y = mvtnorm::rmvnorm(n, mean=rep(0,p), sigma=sig.alpha)
} else if (htype=="H1"){
output$X = mvtnorm::rmvnorm(n, mean=rep(0,p), sigma=sig.alpha)
output$Y = mvtnorm::rmvnorm(n, mean=mu.alpha, sigma=sig.alpha)
}
return(output)
}
# (02) scenario.2013CLX ---------------------------------------------------
#' @keywords internal
#' @noRd
scenario.2013CLX <- function(n, p, no.cov=1, htype=c("H0","H1")){
# covariance modeling 1
scenario.2013CLX.cov1 <- function(p){
Dhalf = diag(sqrt(runif(p, min=0.5, max=2.5)))
SStar = array(0,c(p,p))
for (k in 1:floor(p/5)){
lowk = 5*(k-1)+1
uppk = 5*k
for (i in lowk:uppk){
for (j in lowk:uppk){
SStar[i,j] = 0.5
}
}
}
diag(SStar) = 1.0
mSig = Dhalf%*%SStar%*%Dhalf
return(mSig)
}
scenario.2013CLX.cov2 <- function(p){
Dhalf = diag(sqrt(runif(p, min=0.5, max=2.5)))
sstar = array(0,c(p,p))
for (i in 1:p){
for (j in 1:p){
sstar[i,j] = (0.5^abs(i-j))
}
}
msig = Dhalf%*%sstar%*%Dhalf
return(msig)
}
scenario.2013CLX.cov3 <- function(p){
Dhalf = diag(sqrt(runif(p, min=0.5, max=2.5)))
sstar = array(0,c(p,p))
for (i in 1:(p-1)){
for (j in (i+1):p){
theval = ifelse(rbinom(1, 1, 0.05) > 0.5, 0.5, 0)
sstar[i,j] = theval
sstar[j,i] = theval
}
}
diag(sstar) = 1
del = abs(min(eigen(sstar)$values))+0.05
ssig = (Dhalf%*%(sstar+del*diag(p))%*%Dhalf)/(1+del)
return(ssig)
}
scenario.2013CLX.cov4 <- function(p){
omat = diag(runif(p, min=1, max=5))
tmat = array(0,c(p,p))
for (i in 1:p){
for (j in 1:p){
tmat[i,j] = ((-1)^(i+j))*(0.4^((abs(i-j))^0.1))
}
}
ssig = omat%*%tmat%*%omat
return(ssig)
}
# 1. generate common covariance
no.cov = as.integer(no.cov)
if (!(no.cov %in% 1:4)){
stop("* scenario.2013CLX : 'no.cov' is from 1 to 4.")
}
thesig = switch(no.cov,
"1" = scenario.2013CLX.cov1(p),
"2" = scenario.2013CLX.cov2(p),
"3" = scenario.2013CLX.cov3(p),
"4" = scenario.2013CLX.cov4(p))
# 2. generate samples accordingly
output = list()
htype = match.arg(htype)
if (htype=="H0"){
output$X = mvtnorm::rmvnorm(n, mean=rep(0,p), sigma=thesig)
output$Y = mvtnorm::rmvnorm(n, mean=rep(0,p), sigma=thesig)
} else {
n.upptri = sum(1:(p-1))
vec.utri = rep(0,n.upptri)
vec.utri[sample(1:n.upptri, 4, replace=FALSE)] = max(diag(thesig))*runif(4, min=0, max=4)
Umat = array(0,c(p,p))
Umat[upper.tri(Umat)] = vec.utri
Umat = (Umat + t(Umat))
del = abs(min(c(min(eigen(thesig+Umat)$values), min(eigen(thesig)$values))))+0.05
sig1 = thesig + del*diag(p)
sig2 = thesig + Umat + del*diag(p)
output$X = mvtnorm::rmvnorm(n, mean=rep(0,p), sigma=sig1)
output$Y = mvtnorm::rmvnorm(n, mean=rep(0,p), sigma=sig2)
}
return(output)
}
# (03) gen.cov.2014CLX ----------------------------------------------------
#' @keywords internal
#' @noRd
gen.cov.2014CLX <- function(p, no.cov=1){
scenario.2014CLX.cov1 = function(n, p){
Sig = array(0,c(p,p))
for (k in 1:floor(p/2)){
lowk = as.integer(2*(k-1)+1)
uppk = as.integer(2*k)
for (i in lowk:uppk){
for (j in lowk:uppk){
Sig[i,j] = 0.8
}
}
}
diag(Sig) = 1
return(Sig)
}
scenario.2014CLX.cov2 = function(n, p){
Sig = array(0,c(p,p))
for (i in 1:p){
for (j in 1:p){
Sig[i,j] = 0.6^(abs(i-j))
}
}
return(Sig)
}
scenario.2014CLX.cov3 = function(n, p){
Omega = array(0,c(p,p))
for (i in 1:(p-1)){
Omega[i,i+1] = 0.8
}
for (i in 1:(p-2)){
Omega[i,i+2] = 0.4
}
for (i in 1:(p-3)){
Omega[i,i+3] = 0.2
}
Omega = Omega + t(Omega)
diag(Omega) = 2
Sig = solve(Omega, diag(p))
return(Sig)
}
scenario.2014CLX.cov4 = function(n, p){
Omega = array(0,c(p,p))
for (i in 1:(p-1)){
for (j in (i+1):p){
theval = 0.6^(abs(i-j))
Omega[i,j] = theval
Omega[j,i] = theval
}
}
diag(Omega) = 1
Dhalf = diag(sqrt(runif(p, min=1, max=3)))
Sig = Dhalf%*%solve(Omega,Dhalf)
return(Sig)
}
scenario.2014CLX.cov5 = function(n, p){
Dhalf = diag(sqrt(runif(p, min=1, max=3)))
Ohalf = array(0,c(p,p))
for (k in 1:floor(p/2)){
lowk = as.integer(2*(k-1)+1)
uppk = as.integer(2*k)
for (i in lowk:uppk){
for (j in lowk:uppk){
Ohalf[i,j] = 0.8
}
}
}
diag(Ohalf) = 1
Omega = Dhalf%*%Ohalf%*%Ohalf%*%Dhalf
Sig = solve(Omega, diag(p))
}
scenario.2014CLX.cov6 = function(n, p){
SigStar = array(0,c(p,p))
for (k in 1:floor(p/2)){
lowk = as.integer(2*(k-1)+1)
uppk = as.integer(2*k)
for (i in lowk:uppk){
for (j in lowk:uppk){
SigStar[i,j] = 0.8
}
}
}
diag(SigStar) = 1
Dhalf = diag(sqrt(runif(p, min=1, max=3)))
E = array(0,c(p,p))
noff = sum(1:(p-1))
E[lower.tri(E)] = rbinom(noff, 1, 0.3)*runif(noff, min=-0.2, max=0.2)
E = E + t(E)
triplet = Dhalf%*%SigStar%*%Dhalf
delta = abs(min(eigen(triplet + E)$values))
Sig = triplet + E + delta*diag(p)
return(Sig)
}
scenario.2014CLX.cov7 = function(n, p){
SigStar = array(0,c(p,p))
for (i in 1:(p-1)){
for (j in (i+1):p){
theval = ((abs(i-j))^{-5})/2
SigStar[i,j] = theval
SigStar[j,i] = theval
}
}
diag(SigStar) = 1
Dhalf = diag(sqrt(runif(p, min=1, max=3)))
Sig = Dhalf%*%SigStar%*%Dhalf
return(Sig)
}
scenario.2014CLX.cov8 = function(n, p){
Dhalf = diag(sqrt(runif(p, min=1, max=3)))
Fmat = diag(p)
for (i in 1:(p-1)){
j = i+1
Fmat[i,j] = 0.5
Fmat[j,i] = 0.5
}
UmatOriginal = matrix(rnorm(p*p), nrow=p)
UmatOriginal = UmatOriginal%*%t(UmatOriginal)
U3 = qr.Q(qr(UmatOriginal))[,1:3]
Sig = Dhalf%*%(Fmat + (U3%*%t(U3)))%*%Dhalf
return(Sig)
}
sig.alpha = switch(no.cov,
"1" = scenario.2014CLX.cov1(n,p),
"2" = scenario.2014CLX.cov2(n,p),
"3" = scenario.2014CLX.cov3(n,p),
"4" = scenario.2014CLX.cov4(n,p),
"5" = scenario.2014CLX.cov5(n,p),
"6" = scenario.2014CLX.cov6(n,p),
"7" = scenario.2014CLX.cov7(n,p),
"8" = scenario.2014CLX.cov8(n,p))
return(sig.alpha)
}
# (04) gen.2013CLX.mine ---------------------------------------------------
#' @keywords internal
#' @noRd
gen.2013CLX.mine <- function(n, p, no.cov=1, m, rho){
# m : number of non-zero elements in U
# rho : non-zero signal strength
# covariance modeling 1
scenario.2013CLX.cov1 <- function(p){
Dhalf = diag(sqrt(runif(p, min=0.5, max=2.5)))
SStar = array(0,c(p,p))
for (k in 1:floor(p/5)){
lowk = 5*(k-1)+1
uppk = 5*k
for (i in lowk:uppk){
for (j in lowk:uppk){
SStar[i,j] = 0.5
}
}
}
diag(SStar) = 1.0
mSig = Dhalf%*%SStar%*%Dhalf
return(mSig)
}
scenario.2013CLX.cov2 <- function(p){
Dhalf = diag(sqrt(runif(p, min=0.5, max=2.5)))
sstar = array(0,c(p,p))
for (i in 1:p){
for (j in 1:p){
sstar[i,j] = (0.5^abs(i-j))
}
}
msig = Dhalf%*%sstar%*%Dhalf
return(msig)
}
scenario.2013CLX.cov3 <- function(p){
Dhalf = diag(sqrt(runif(p, min=0.5, max=2.5)))
sstar = array(0,c(p,p))
for (i in 1:(p-1)){
for (j in (i+1):p){
theval = ifelse(rbinom(1, 1, 0.05) > 0.5, 0.5, 0)
sstar[i,j] = theval
sstar[j,i] = theval
}
}
diag(sstar) = 1
del = abs(min(eigen(sstar)$values))+0.05
ssig = (Dhalf%*%(sstar+del*diag(p))%*%Dhalf)/(1+del)
return(ssig)
}
scenario.2013CLX.cov4 <- function(p){
omat = diag(runif(p, min=1, max=5))
tmat = array(0,c(p,p))
for (i in 1:p){
for (j in 1:p){
tmat[i,j] = ((-1)^(i+j))*(0.4^((abs(i-j))^0.1))
}
}
ssig = omat%*%tmat%*%omat
return(ssig)
}
# 1. generate common covariance
no.cov = as.integer(no.cov)
if (!(no.cov %in% 1:4)){
stop("* scenario.2013CLX : 'no.cov' is from 1 to 4.")
}
thesig = switch(no.cov,
"1" = scenario.2013CLX.cov1(p),
"2" = scenario.2013CLX.cov2(p),
"3" = scenario.2013CLX.cov3(p),
"4" = scenario.2013CLX.cov4(p))
# 2. generate samples accordingly
output = list()
# 2-1. Let's adjust for the signal strength
n.upptri = sum(1:(p-1))
vec.utri = rep(0,n.upptri)
vec.utri[sample(1:n.upptri, m, replace=FALSE)] = rho # fixed strength
Umat = array(0,c(p,p))
Umat[upper.tri(Umat)] = vec.utri
Umat = (Umat + t(Umat))
del = abs(min(c(min(eigen(thesig+Umat)$values), min(eigen(thesig)$values))))+0.05
sig1 = thesig + del*diag(p)
sig2 = thesig + Umat + del*diag(p)
# 2-2. generate data for H0
tmpH0 = list()
tmpH0$X = mvtnorm::rmvnorm(n, mean=rep(0,p), sigma=sig1)
tmpH0$Y = mvtnorm::rmvnorm(n, mean=rep(0,p), sigma=sig1)
# 2-3. generate data for H1
tmpH1 = list()
tmpH1$X = mvtnorm::rmvnorm(n, mean=rep(0,p), sigma=sig1)
tmpH1$Y = mvtnorm::rmvnorm(n, mean=rep(0,p), sigma=sig2)
# 2-4. wrap-up and return
output$H0 = tmpH0
output$H1 = tmpH1
return(output)
}
# (05) gen.2013CLX.dense --------------------------------------------------
#' @keywords internal
#' @noRd
gen.2013CLX.dense <- function(n, p, no.cov=1, rho){
# rho : non-zero signal strength
# covariance modeling 1
scenario.2013CLX.cov1 <- function(p){
Dhalf = diag(sqrt(runif(p, min=0.5, max=2.5)))
SStar = array(0,c(p,p))
for (k in 1:floor(p/5)){
lowk = 5*(k-1)+1
uppk = 5*k
for (i in lowk:uppk){
for (j in lowk:uppk){
SStar[i,j] = 0.5
}
}
}
diag(SStar) = 1.0
mSig = Dhalf%*%SStar%*%Dhalf
return(mSig)
}
scenario.2013CLX.cov2 <- function(p){
Dhalf = diag(sqrt(runif(p, min=0.5, max=2.5)))
sstar = array(0,c(p,p))
for (i in 1:p){
for (j in 1:p){
sstar[i,j] = (0.5^abs(i-j))
}
}
msig = Dhalf%*%sstar%*%Dhalf
return(msig)
}
scenario.2013CLX.cov3 <- function(p){
Dhalf = diag(sqrt(runif(p, min=0.5, max=2.5)))
sstar = array(0,c(p,p))
for (i in 1:(p-1)){
for (j in (i+1):p){
theval = ifelse(rbinom(1, 1, 0.05) > 0.5, 0.5, 0)
sstar[i,j] = theval
sstar[j,i] = theval
}
}
diag(sstar) = 1
del = abs(min(eigen(sstar)$values))+0.05
ssig = (Dhalf%*%(sstar+del*diag(p))%*%Dhalf)/(1+del)
return(ssig)
}
scenario.2013CLX.cov4 <- function(p){
omat = diag(runif(p, min=1, max=5))
tmat = array(0,c(p,p))
for (i in 1:p){
for (j in 1:p){
tmat[i,j] = ((-1)^(i+j))*(0.4^((abs(i-j))^0.1))
}
}
ssig = omat%*%tmat%*%omat
return(ssig)
}
# 1. generate common covariance
no.cov = as.integer(no.cov)
if (!(no.cov %in% 1:4)){
stop("* scenario.2013CLX : 'no.cov' is from 1 to 4.")
}
thesig = switch(no.cov,
"1" = scenario.2013CLX.cov1(p),
"2" = scenario.2013CLX.cov2(p),
"3" = scenario.2013CLX.cov3(p),
"4" = scenario.2013CLX.cov4(p))
# 2. generate samples accordingly
output = list()
# 2-1. Let's adjust for the signal strength by Kyoungjae's method
vvec = rnorm(p)
vvec = vvec/sqrt(sum(vvec^2))
sig1 = thesig
sig2 = thesig + outer(vvec,vvec)*rho
# 2-2. generate data for H0
tmpH0 = list()
tmpH0$X = mvtnorm::rmvnorm(n, mean=rep(0,p), sigma=sig1)
tmpH0$Y = mvtnorm::rmvnorm(n, mean=rep(0,p), sigma=sig1)
# 2-3. generate data for H1
tmpH1 = list()
tmpH1$X = mvtnorm::rmvnorm(n, mean=rep(0,p), sigma=sig1)
tmpH1$Y = mvtnorm::rmvnorm(n, mean=rep(0,p), sigma=sig2)
# 2-4. wrap-up and return
output$H0 = tmpH0
output$H1 = tmpH1
return(output)
}
#' @keywords internal
#' @noRd
gen.2013CLX.dense.ones <- function(n, p, no.cov=1, rho){
# rho : non-zero signal strength
# covariance modeling 1
scenario.2013CLX.cov1 <- function(p){
Dhalf = diag(sqrt(runif(p, min=0.5, max=2.5)))
SStar = array(0,c(p,p))
for (k in 1:floor(p/5)){
lowk = 5*(k-1)+1
uppk = 5*k
for (i in lowk:uppk){
for (j in lowk:uppk){
SStar[i,j] = 0.5
}
}
}
diag(SStar) = 1.0
mSig = Dhalf%*%SStar%*%Dhalf
return(mSig)
}
scenario.2013CLX.cov2 <- function(p){
Dhalf = diag(sqrt(runif(p, min=0.5, max=2.5)))
sstar = array(0,c(p,p))
for (i in 1:p){
for (j in 1:p){
sstar[i,j] = (0.5^abs(i-j))
}
}
msig = Dhalf%*%sstar%*%Dhalf
return(msig)
}
scenario.2013CLX.cov3 <- function(p){
Dhalf = diag(sqrt(runif(p, min=0.5, max=2.5)))
sstar = array(0,c(p,p))
for (i in 1:(p-1)){
for (j in (i+1):p){
theval = ifelse(rbinom(1, 1, 0.05) > 0.5, 0.5, 0)
sstar[i,j] = theval
sstar[j,i] = theval
}
}
diag(sstar) = 1
del = abs(min(eigen(sstar)$values))+0.05
ssig = (Dhalf%*%(sstar+del*diag(p))%*%Dhalf)/(1+del)
return(ssig)
}
scenario.2013CLX.cov4 <- function(p){
omat = diag(runif(p, min=1, max=5))
tmat = array(0,c(p,p))
for (i in 1:p){
for (j in 1:p){
tmat[i,j] = ((-1)^(i+j))*(0.4^((abs(i-j))^0.1))
}
}
ssig = omat%*%tmat%*%omat
return(ssig)
}
# 1. generate common covariance
no.cov = as.integer(no.cov)
if (!(no.cov %in% 1:4)){
stop("* scenario.2013CLX : 'no.cov' is from 1 to 4.")
}
thesig = switch(no.cov,
"1" = scenario.2013CLX.cov1(p),
"2" = scenario.2013CLX.cov2(p),
"3" = scenario.2013CLX.cov3(p),
"4" = scenario.2013CLX.cov4(p))
# 2. generate samples accordingly
output = list()
# 2-1. Let's adjust for the signal strength by Kyoungjae's method
vvec = rep(1,p)
sig1 = thesig
sig2 = thesig + outer(vvec,vvec)*rho
# 2-2. generate data for H0
tmpH0 = list()
tmpH0$X = mvtnorm::rmvnorm(n, mean=rep(0,p), sigma=sig1)
tmpH0$Y = mvtnorm::rmvnorm(n, mean=rep(0,p), sigma=sig1)
# 2-3. generate data for H1
tmpH1 = list()
tmpH1$X = mvtnorm::rmvnorm(n, mean=rep(0,p), sigma=sig1)
tmpH1$Y = mvtnorm::rmvnorm(n, mean=rep(0,p), sigma=sig2)
# 2-4. wrap-up and return
output$H0 = tmpH0
output$H1 = tmpH1
return(output)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.