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)
|