# m is the rank of supspace (generically should be group size) # Z is the observed norm # finalinterval is R_Y # alpha is the level controled # maxL is the maximum considered L. # output: L (left endpoint of CI) BinarySearch=function(m,Z,finalintervals,alpha,tol=0.01,maxL=2*Z){ err=1+tol; L_Upper=maxL; L_Lower=-m; fl=length(finalintervals); leftend=finalintervals[seq(1,fl,2)]; rightend=finalintervals[seq(2,fl,2)]; ll=ceiling(fl/2); f1=rep(0,ll); f2=rep(0,ll); temp=which(leftend<=Z); posi=temp[length(temp)]; while (err>tol){ L=(L_Upper+L_Lower)/2; integrand<-function(r) {r^(m-1)*exp(-(r-L)^2/2)} for (i in 1:ll){ f1[i]=integrate(integrand,leftend[i],rightend[i])\$value; if (i=0) || is.nan(f) ){ L_Upper=L; } else { L_Lower=L; } err=L_Upper-L_Lower; } return(L) } # Construct union of finite intervals within [lower,upper] such that points # between lowerset(i) and upperset(i) are excluded for all i. unioninterval=function(lowerset,upperset,lower,upper){ ll=length(lowerset); finterval=c(lower,upper); for (i in 1:ll){ left=lowerset[i]; right=upperset[i]; temp1=which(finterval>=left); temp2=which(finterval<=right); if ((length(temp1)==0) || (length(temp2)==0)){ next } index1=temp1[1]; index2=temp2[length(temp2)]; if (index1>1){ LL=finterval[1:(index1-1)]; } else { LL=c(); } if (index2=||c2*Z+d2|| a=sum(c1^2)-sum(c2^2); b=2*(sum(c1*d1)-sum(c2*d2)); c=sum(d1^2)-sum(d2^2); if (a<0){ #print(b^2-4*a*c) quadroot_plus = (-b+sqrt(b^2-4*a*c))/(2*a) quadroot_minus = (-b-sqrt(b^2-4*a*c))/(2*a) Z_lower = max(Z_lower,quadroot_plus) Z_upper = min(Z_upper,quadroot_minus) } if (a>0){ if (b^2-4*a*c>0){ quadroot_plus = (-b+sqrt(b^2-4*a*c))/(2*a); quadroot_minus = (-b-sqrt(b^2-4*a*c))/(2*a); upperset=c(upperset,quadroot_plus); lowerset=c(lowerset,quadroot_minus); } } if(a==0){ if (b>0){ Z_lower=max(Z_lower,-c/b); } if (b<-0){ Z_upper=min(Z_upper,-c/b); } } } } # updata ct and dt ct[setdiff(1:(G*m),position)]=0; dt[setdiff(1:(G*m),position)]=0 ct=-(lr/n)*t(X)%*%u+(lr/n*XtX+diag(G*m))%*%ct; dt=-(lr/n)*t(X)%*%yperp+(lr/n*XtX+diag(G*m))%*%dt; } finalintervals=unioninterval(lowerset,upperset,Z_lower,Z_upper)/sigma; Z=Z/sigma; # do one-sided p-value fl=length(finalintervals); leftend=finalintervals[seq(1,fl,2)]; rightend=finalintervals[seq(2,fl,2)]; denominator=sum(pchisq(rightend^2,m)-pchisq(leftend^2,m)); if (denominator==0){ pval[g0]=0; } else { temp=which(leftend<=Z); posi=temp[length(temp)]; numerator=sum(pchisq(rightend[posi:length(rightend)]^2,m)-pchisq(leftend[posi:length(leftend)]^2,m)) - (pchisq(Z^2,m)-pchisq(leftend[posi]^2,m)); p1=numerator/denominator; pval[g0]=p1; } # do one-sided confidence interval L=sigma*BinarySearch(m,Z,finalintervals,alpha); CI[g0]=L; } return(list('pval'=pval,'CI'=CI)) }