Author:

Zhong, P-S, Chen, S. X. and Xu M.

Title:

Zhong, P-S, Chen, S. X. and Xu M. (2013). Tests alternative to higher criticism for high dimensional means under sparsity and column-wise dependence, Annals of Statistics, 41, 2820-2851.
Gaussvec<-function(m, p, n, rho, mu)
{
     x=matrix(0,p,n)
     lam=rep(0,m)
     lam[1:(m/2+1)]=rho^(0:(m/2))
     lam[(m/2+2):m]=rho^((m/2-1):1)
     lam0=fft(lam);
     
     urv=matrix(rnorm(m*n,0,1),m,n)
     aseq=matrix(0,m,n)
     aseq[1,]=sqrt(lam0[1])*urv[1,]/sqrt(m);
     aseq[m/2+1,]=sqrt(lam0[m/2+1])*urv[m/2+1,]/sqrt(m);
     aseq[2:(m/2),]=sqrt(lam0[2:(m/2)])*(urv[2:(m/2),]+1i*urv[(2+m/2):m,])/sqrt(2*m);
     aseq[m:(m/2+2),]=sqrt(lam0[2:(m/2)])*(urv[2:(m/2),]-1i*urv[(2+m/2):m,])/sqrt(2*m)
     aseq0=mvfft(aseq);
     x=Re(aseq0)[1:p,]+mu    
     return(x);
}

qvalue<-function(pvalues, p)
{
   pvalues=p*sort(pvalues)/c(1:p);
   for (s in 1:p)
   {pvalues[s]=min(pvalues[s:p])}
   return(pvalues);
}


TestStats<-function(X,eta)
{
p<-dim(X)[1]
n<-dim(X)[2]
Xbar<-apply(X,1,mean)

tstat<-abs(sqrt(n)*Xbar) 
pvalues<-2*(1-pnorm(tstat))
FDRpvalues<-qvalue(pvalues,p)[1]

lambda2<-sort(n*(Xbar)^2)
lamind<-(lambda2<2*(1-eta)*log(p))*c(1:p)
lambda<-sqrt(lambda2)    
muT2<-p*(2*lambda[lamind]*dnorm(lambda[lamind])+2*pnorm(lambda[lamind],0,1,lower.tail=FALSE))
varT2<-p*(2*(lambda[lamind]^3+3*lambda[lamind])*dnorm(lambda[lamind])+6*pnorm(lambda[lamind],0,1,lower.tail=FALSE))-(muT2^2)/p
Ynvec<-n*(Xbar)^2
sortYn<-sort(Ynvec,decreasing=TRUE)
Tn2<-sort(cumsum(sortYn),decreasing=TRUE)
stdTn2<-(Tn2[lamind]-muT2)/sqrt(varT2)
maxstdTn2<-max(stdTn2)

muT1<-p*2*dnorm(lambda[lamind])
varT1<-p*(2*lambda[lamind]*dnorm(lambda[lamind])+2*pnorm(lambda[lamind],0,1,lower.tail=FALSE))-(muT1^2)/p
Tn1<-sort(cumsum(sqrt(sortYn)),decreasing=TRUE)
stdTn1<-(Tn1[lamind]-muT1)/sqrt(varT1)
maxstdTn1<-max(stdTn1)

muHC<-2*p*pnorm(lambda[lamind],0,1,lower.tail=FALSE)
varHC<-2*p*pnorm(lambda[lamind],0,1,lower.tail=FALSE)*(1-2*pnorm(lambda[lamind],0,1,lower.tail=FALSE))
cumHC<-sort(cumsum(rep(1,p)),decreasing=TRUE)
HC<-(cumHC[lamind]-muHC)/sqrt(varHC)
maxHC<-max(HC)

return(list(HC=maxHC,L1=maxstdTn1,L2=maxstdTn2,FDR=FDRpvalues))
}



powcomp<-function(p,n,beta,r,rho,reps)
{
g<-trunc((log(p)/log(2))+1)+1
m<-2^g
nonzeros<-floor(p^(1-beta))
mus<-c(rep(sqrt(2*r*log(p)/n),nonzeros),rep(0,p-nonzeros))
HCvec<-maxTnvec2<-maxTnvec1<-HCstdvec<-maxTnstdvec2<-maxTnstdvec1<-FDRpvalues<-NULL
eta<-0.05
an<-sqrt(2*log(log(p)))
bn<-2*log(log(p))+0.5*log(log(log(p)))-0.5*log(4*pi/((1-eta)^2))

for (j in 1:reps)
{
 X<-Gaussvec(m, p, n, rho, mus)
 results<-TestStats(X,0.05)
 HCvec[j]<-results$HC
 maxTnvec1[j]<-results$L1
 maxTnvec2[j]<-results$L2
 FDRpvalues[j]<-results$FDR
}
 HCstdvec<-an*HCvec-bn
 maxTnstdvec1<-an*maxTnvec1-bn
 maxTnstdvec2<-an*maxTnvec2-bn
return(list(HC=HCvec,L1=maxTnvec1,L2=maxTnvec2,stdHC=HCstdvec,stdL1=maxTnstdvec1,stdL2=maxTnstdvec2,FDR=FDRpvalues))
}

exq<-function(p){return(-log(-log(p)))}

## The following function compute the test statistics proposed in ZCX (2013, AoS)
## and their comparison with Higher Criticism (Donoho and Jin, 2004) and the FDR based method.
## It was used in simulation setting, which can be modified for your purpose.
## p: data dimension
## n: sample size
## beta: sparsity parameter
## r: signal level
## rho: correlation in Gaussian random vector
## reps: the number of simulation replications

## One example is given below

nullp2500n40rho5<-powcomp(2500,40,0.7,0,0.5,50)