#############################################################################################


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