############################################################################################# # This script computes k-dimensional Quasi-Conventional simultaneous confidence intervals # # for a specified sample point x=(x_1,x_2,...,x_k) # # This is a multivariate generalization # to the Univarate Quasi-Conventional intervals # with more power to determine sign # Benjamini, Hochberg and Stark (1998) # # ###################################################################################### # Parameters: # # r is a value larger than 1 that controls the maximal length<=2rcalpha # # clevel is the significance level (default=0.95) # # x is the desired sample vector (for example x=c(3,4)) # # # ###################################################################################### ###Examples #QCCI(x=c(1.9473,2.1335),r=1.2,clevel=.95) #QCCI(x=c(1.9473,2.1335,2.558),r=1.8,clevel=.95) ###################################################################################### # Warning !!! this is perliminal script # its validity needs to be rechecked. ################################################################ QCCI<-function(x,r=1.2,clevel=.95){ dim<-length(x) xord<-order(abs(x)) sx<-sign(x) calpha<-qnorm(0.5*(exp(log(clevel)/dim)+1)) rcalpha<-r*calpha chi<-rep(0,dim) for (i in (1:dim)) { rclevel<-clevel/(2*pnorm(rcalpha)-1)**(dim-i) rclevel<-exp(log(rclevel)/i) c<-calpha for (itr in 1:10) { u<-(pnorm(2*rcalpha-c)+pnorm(c)-1 )-rclevel v<-dnorm(c)-dnorm(2*rcalpha-c) diff<-abs(u/v) c<-c-u/v if (diff<0.000001) break } chi[i]<-min(c,2*rcalpha-c) } QC_low <-abs(x)-rcalpha QC_up <-abs(x)+rcalpha if (min(abs(x))>2*rcalpha) { QC_low <-abs(x)-calph QC_up <-abs(x)+calpha } for (l in (1:dim)) { for (a in (1:dim)) { if (xord[a]<=l && abs(x[a])<=2*rcalpha) { QC_low[a]<-max(0,abs(x[a])-2*rcalpha+chi[1]) QC_low[a]<-QC_low[a]+min(abs(x[a])-chi[l],0) if (max(abs(x))>2*rcalpha) { QC_up[a]<-abs(x[a])+calpha } }}} if(max(abs(x))<=chi[1]){QC_low<- x-rcalpha} QCL <- QC_low*sx QCU <- QC_up*sx QCCI<-rbind(x,QCL,QCU) return(QCCI,chi) } ###Examples QCCI(x=c(1.9473,2.1335),r=1.2,clevel=.95) QCCI(x=c(1.947,2.134,2.558),r=1.2,clevel=.95) QCCI(x=c(1.947,2.134,-2.060,2.558),r=2,clevel=.95)