### R.code for SCIs for all diff CVs in DLN (k=3)
library(EnvStats)
library(HDInterval)
library(invgamma)
library(plotrix)
CI.SCIdiffCVs <- function(M=100,m=5000,n1=25,n2=25,n3=25,var1=0.5,var2=0.5,var3=0.5,mean1=0,mean2=0,mean3=0,
                     delta1=0.5,delta2=0.5,delta3=0.5,alpha=0.05){
starttime <- proc.time()
  
  cv1 <- sqrt(exp(var1)-1)                                  #variance of the LN distribution in population 1
  cv2 <- sqrt(exp(var2)-1)                                  #variance of the LN distribution in population 2
  cv3 <- sqrt(exp(var3)-1)                                  #variance of the LN distribution in population 3
  ex1 <- delta1*exp(mean1+var1/2)                           #mean of the delta-LN distribution in population 1
  ex2 <- delta2*exp(mean2+var2/2)                           #mean of the delta-LN distribution in population 2
  ex3 <- delta3*exp(mean3+var3/2)                           #mean of the delta-LN distribution in population 3
  var_x1 <- delta1*exp((2*mean1)+var1)*(exp(var1)-delta1)   #variance of the delta-LN distribution in population 1
  var_x2 <- delta2*exp((2*mean2)+var2)*(exp(var2)-delta2)   #variance of the delta-LN distribution in population 2
  var_x3 <- delta3*exp((2*mean3)+var3)*(exp(var3)-delta3)   #variance of the delta-LN distribution in population 3
  eta1 <- sqrt(var_x1)/ex1                                  #coefficient of variation of the delta-LN distribution in population 1
  eta2 <- sqrt(var_x2)/ex2                                  #coefficient of variation of the delta-LN distribution in population 2
  eta3 <- sqrt(var_x3)/ex3                                  #coefficient of variation of the delta-LN distribution in population 3
  gamma12 <- eta1-eta2                                      #difference between eta1 and eta2
  gamma13 <- eta1-eta3                                      #difference between eta1 and eta3
  gamma23 <- eta2-eta3                                      #difference between eta2 and eta3
  
  #FGCI
  U1 <- matrix(rep(0,m))
  U2 <- matrix(rep(0,m))
  U3 <- matrix(rep(0,m))
  beta11 <- matrix(rep(0,m))
  beta12 <- matrix(rep(0,m))
  beta21 <- matrix(rep(0,m))
  beta22 <- matrix(rep(0,m))
  beta31 <- matrix(rep(0,m))
  beta32 <- matrix(rep(0,m))
  var1_GFQ <- matrix(rep(0,m))
  var2_GFQ <- matrix(rep(0,m))
  var3_GFQ <- matrix(rep(0,m))
  delta1_GFQ <- matrix(rep(0,m))
  delta2_GFQ <- matrix(rep(0,m))
  delta3_GFQ <- matrix(rep(0,m))
  R1.GFQ <- matrix(rep(0,m))
  R2.GFQ <- matrix(rep(0,m))
  R3.GFQ <- matrix(rep(0,m))
  R12.GFQ <- matrix(rep(0,m))
  R13.GFQ <- matrix(rep(0,m))
  R23.GFQ <- matrix(rep(0,m))
  
  #Bayesian 
  var1.jrule <- matrix(rep(0,m))
  var2.jrule <- matrix(rep(0,m))
  var3.jrule <- matrix(rep(0,m))
  delta1.jrule <- matrix(rep(0,m))
  delta2.jrule <- matrix(rep(0,m))
  delta3.jrule <- matrix(rep(0,m))
  beta1.jrule <- matrix(rep(0,m))
  beta2.jrule <- matrix(rep(0,m))
  beta3.jrule <- matrix(rep(0,m))
  eta1.jrule <- matrix(rep(0,m))
  eta2.jrule <- matrix(rep(0,m))
  eta3.jrule <- matrix(rep(0,m))
  gamma12.jrule <- matrix(rep(0,m))
  gamma13.jrule <- matrix(rep(0,m))
  gamma23.jrule <- matrix(rep(0,m))
  
  var1.uni <- matrix(rep(0,m))
  var2.uni <- matrix(rep(0,m))
  var3.uni <- matrix(rep(0,m))
  delta1.uni <- matrix(rep(0,m))
  delta2.uni <- matrix(rep(0,m))
  delta3.uni <- matrix(rep(0,m))
  beta1.uni <- matrix(rep(0,m))
  beta2.uni <- matrix(rep(0,m))
  beta3.uni <- matrix(rep(0,m))
  eta1.uni <- matrix(rep(0,m))
  eta2.uni <- matrix(rep(0,m))
  eta3.uni <- matrix(rep(0,m))
  gamma12.uni <- matrix(rep(0,m))
  gamma13.uni <- matrix(rep(0,m))
  gamma23.uni <- matrix(rep(0,m))
  
  #FGCI
  L.rFGCI12 <- numeric()
  U.rFGCI12 <- numeric()
  CP.rFGCI12 <- numeric()
  Length.rFGCI12 <- numeric()
  L.rFGCI13 <- numeric()
  U.rFGCI13 <- numeric()
  CP.rFGCI13 <- numeric()
  Length.rFGCI13 <- numeric()
  L.rFGCI23 <- numeric()
  U.rFGCI23 <- numeric()
  CP.rFGCI23 <- numeric()
  Length.rFGCI23 <- numeric()
  CP.rFGCI <- numeric()
  AL.rFGCI <- numeric()
  SE.Length.rFGCI <- numeric()
  
  #Equitailed Bayesian based on Jeffreys' rule prior
  L.rBaye12_jrule <- numeric()
  U.rBaye12_jrule <- numeric()
  CP.rBaye12_jrule <- numeric()
  Length.rBaye12_jrule <- numeric()
  L.rBaye13_jrule <- numeric()
  U.rBaye13_jrule <- numeric()
  CP.rBaye13_jrule <- numeric()
  Length.rBaye13_jrule <- numeric()
  L.rBaye23_jrule <- numeric()
  U.rBaye23_jrule <- numeric()
  CP.rBaye23_jrule <- numeric()
  Length.rBaye23_jrule <- numeric()
  CP.rBaye_jrule <- numeric()
  AL.rBaye_jrule <- numeric()
  SE.Length.rBaye_jrule <- numeric()
  
  #Equitailed Bayesian based on uniform prior
  L.rBaye12_uni <- numeric()
  U.rBaye12_uni <- numeric()
  CP.rBaye12_uni <- numeric()
  Length.rBaye12_uni <- numeric()
  L.rBaye13_uni <- numeric()
  U.rBaye13_uni <- numeric()
  CP.rBaye13_uni <- numeric()
  Length.rBaye13_uni <- numeric()
  L.rBaye23_uni <- numeric()
  U.rBaye23_uni <- numeric()
  CP.rBaye23_uni <- numeric()
  Length.rBaye23_uni <- numeric()
  CP.rBaye_uni <- numeric()
  AL.rBaye_uni <- numeric()
  SE.Length.rBaye_uni <- numeric() 
  
  #HPD Bayesian based on Jeffreys' rule prior
  hpd12.jrule <- matrix(rep(0,2))
  L.rBaye12_hpd.jrule <- numeric() 
  U.rBaye12_hpd.jrule <- numeric() 
  CP.rBaye12_hpd.jrule <- numeric() 
  Length.rBaye12_hpd.jrule <- numeric() 
  hpd13.jrule <- matrix(rep(0,2))
  L.rBaye13_hpd.jrule <- numeric() 
  U.rBaye13_hpd.jrule <- numeric() 
  CP.rBaye13_hpd.jrule <- numeric() 
  Length.rBaye13_hpd.jrule <- numeric() 
  hpd23.jrule <- matrix(rep(0,2))
  L.rBaye23_hpd.jrule <- numeric() 
  U.rBaye23_hpd.jrule <- numeric() 
  CP.rBaye23_hpd.jrule <- numeric() 
  Length.rBaye23_hpd.jrule <- numeric() 
  CP.rBaye_hpd.jrule <- numeric() 
  AL.rBaye_hpd.jrule <- numeric() 
  SE.Length.rBaye_hpd.jrule <- numeric()
  
  #HPD Bayesian based on uniform prior
  hpd12.uni <- matrix(rep(0,2))
  L.rBaye12_hpd.uni <- numeric() 
  U.rBaye12_hpd.uni <- numeric() 
  CP.rBaye12_hpd.uni <- numeric() 
  Length.rBaye12_hpd.uni <- numeric() 
  hpd13.uni <- matrix(rep(0,2))
  L.rBaye13_hpd.uni <- numeric() 
  U.rBaye13_hpd.uni <- numeric() 
  CP.rBaye13_hpd.uni <- numeric() 
  Length.rBaye13_hpd.uni <- numeric() 
  hpd23.uni <- matrix(rep(0,2))
  L.rBaye23_hpd.uni <- numeric() 
  U.rBaye23_hpd.uni <- numeric() 
  CP.rBaye23_hpd.uni <- numeric() 
  Length.rBaye23_hpd.uni <- numeric() 
  CP.rBaye_hpd.uni <- numeric() 
  AL.rBaye_hpd.uni <- numeric() 
  SE.Length.rBaye_hpd.uni <- numeric()
  
  #MOVER based on Wilson score method
  L.rMOVER12_W <- numeric()
  U.rMOVER12_W <- numeric()
  CP.rMOVER12_W <- numeric()
  Length.rMOVER12_W <- numeric()
  L.rMOVER13_W <- numeric()
  U.rMOVER13_W <- numeric()
  CP.rMOVER13_W <- numeric()
  Length.rMOVER13_W <- numeric()
  L.rMOVER23_W <- numeric()
  U.rMOVER23_W <- numeric()
  CP.rMOVER23_W <- numeric()
  Length.rMOVER23_W <- numeric()
  CP.rMOVER_W <- numeric()
  AL.rMOVER_W <- numeric()
  SE.Length.rMOVER_W <- numeric()
  
  i = 1
  while(i <= M){                                               #Begin loop i
    
    x1 <- rzmlnormAlt(n1,exp(mean1+var1/2),cv1,1-delta1)  
    x2 <- rzmlnormAlt(n2,exp(mean2+var2/2),cv2,1-delta2)
    x3 <- rzmlnormAlt(n3,exp(mean3+var3/2),cv3,1-delta3)
    
    nonzero_value1 <- log(x1[which(x1!=0)])                    
    nonzero_value2 <- log(x2[which(x2!=0)])                    
    nonzero_value3 <- log(x3[which(x3!=0)])                    
    
    size_nonzero1 <- length(nonzero_value1)                    
    if(size_nonzero1 < n1*delta1){
      next 
    }
    
    size_nonzero2 <- length(nonzero_value2)                    
    if(size_nonzero2 < n2*delta2){
      next 
    }
    
    size_nonzero3 <- length(nonzero_value3)                    
    if(size_nonzero3 < n3*delta3){
      next 
    }
    
    mean_hat1 <- mean(nonzero_value1)
    var_hat1 <- var(nonzero_value1)
    delta_hat1 <- size_nonzero1/n1
    eta_hat1 <- sqrt((exp(var_hat1)-delta_hat1)/delta_hat1)
    
    mean_hat2 <- mean(nonzero_value2)
    var_hat2 <- var(nonzero_value2)
    delta_hat2 <- size_nonzero2/n2
    eta_hat2 <- sqrt((exp(var_hat2)-delta_hat2)/delta_hat2)
    
    mean_hat3 <- mean(nonzero_value3)
    var_hat3 <- var(nonzero_value3)
    delta_hat3 <- size_nonzero3/n3
    eta_hat3 <- sqrt((exp(var_hat3)-delta_hat3)/delta_hat3)
    
    #Fiducial generalized confidence interval..................................................................
    U1 <- rchisq(m,(size_nonzero1-1))                                   
    U2 <- rchisq(m,(size_nonzero2-1))                                   
    U3 <- rchisq(m,(size_nonzero3-1))                                   
    beta11 <- rbeta(m,size_nonzero1,(n1-size_nonzero1)+1)                 
    beta12 <- rbeta(m,size_nonzero1+1,n1-size_nonzero1)
    beta21 <- rbeta(m,size_nonzero2,(n2-size_nonzero2)+1)                 
    beta22 <- rbeta(m,size_nonzero2+1,n2-size_nonzero2)
    beta31 <- rbeta(m,size_nonzero3,(n3-size_nonzero3)+1)                 
    beta32 <- rbeta(m,size_nonzero3+1,n3-size_nonzero3)
    var1_GFQ <- ((size_nonzero1-1)*var_hat1)/U1                                       
    var2_GFQ <- ((size_nonzero2-1)*var_hat2)/U2                                       
    var3_GFQ <- ((size_nonzero3-1)*var_hat3)/U3                                       
    delta1_GFQ <- (1/2)*beta11+(1/2)*beta12                                           
    delta2_GFQ <- (1/2)*beta21+(1/2)*beta22                                           
    delta3_GFQ <- (1/2)*beta31+(1/2)*beta32                                           
    R1.GFQ <- sqrt((exp(var1_GFQ)-delta1_GFQ)/delta1_GFQ)                             
    R2.GFQ <- sqrt((exp(var2_GFQ)-delta2_GFQ)/delta2_GFQ)                             
    R3.GFQ <- sqrt((exp(var3_GFQ)-delta3_GFQ)/delta3_GFQ)                             
    R12.GFQ <- R1.GFQ-R2.GFQ                                                          
    if(is.na(R12.GFQ)){next}
    R13.GFQ <- R1.GFQ-R3.GFQ                                                          
    if(is.na(R13.GFQ)){next}
    R23.GFQ <- R2.GFQ-R3.GFQ                                                          
    if(is.na(R23.GFQ)){next}
    
    L.rFGCI12[i] <-quantile(R12.GFQ,alpha/2)                                          
    U.rFGCI12[i] <-quantile(R12.GFQ,(1-(alpha/2)))                                    
    CP.rFGCI12[i] <- ifelse(L.rFGCI12[i]< gamma12 && gamma12 <U.rFGCI12[i],1,0)       
    Length.rFGCI12[i] <- U.rFGCI12[i]-L.rFGCI12[i]                                    
    L.rFGCI13[i] <-quantile(R13.GFQ,alpha/2)                                          
    U.rFGCI13[i] <-quantile(R13.GFQ,(1-(alpha/2)))                                    
    CP.rFGCI13[i] <- ifelse(L.rFGCI13[i]< gamma13 && gamma13 <U.rFGCI13[i],1,0)       
    Length.rFGCI13[i] <- U.rFGCI13[i]-L.rFGCI13[i]                                    
    L.rFGCI23[i] <-quantile(R23.GFQ,alpha/2)                                          
    U.rFGCI23[i] <-quantile(R23.GFQ,(1-(alpha/2)))                                    
    CP.rFGCI23[i] <- ifelse(L.rFGCI23[i]< gamma23 && gamma23 <U.rFGCI23[i],1,0)       
    Length.rFGCI23[i] <- U.rFGCI23[i]-L.rFGCI23[i]                                    
    CP.rFGCI[i] <- (CP.rFGCI12[i]+CP.rFGCI13[i]+CP.rFGCI23[i])/3                      
    AL.rFGCI[i] <- (Length.rFGCI12[i]+Length.rFGCI13[i]+Length.rFGCI23[i])/3                        
    
    #Equitailed Bayesian based on Jeffreys' rule prior.........................................................
    beta1.jrule <- rbeta(m,size_nonzero1+(3/2),(n1-size_nonzero1)+(1/2))                  
    beta2.jrule <- rbeta(m,size_nonzero2+(3/2),(n2-size_nonzero2)+(1/2))                  
    beta3.jrule <- rbeta(m,size_nonzero3+(3/2),(n3-size_nonzero3)+(1/2))                  
    var1.jrule <- rinvgamma(m,size_nonzero1/2,(size_nonzero1*var_hat1)/2,scale = 1/rate)
    var2.jrule <- rinvgamma(m,size_nonzero2/2,(size_nonzero2*var_hat2)/2,scale = 1/rate)
    var3.jrule <- rinvgamma(m,size_nonzero3/2,(size_nonzero3*var_hat3)/2,scale = 1/rate)
    delta1.jrule <- beta1.jrule                                                           
    delta2.jrule <- beta2.jrule                                                           
    delta3.jrule <- beta3.jrule                                                           
    eta1.jrule <- sqrt((exp(var1.jrule)-delta1.jrule)/delta1.jrule)                       
    eta2.jrule <- sqrt((exp(var2.jrule)-delta2.jrule)/delta2.jrule)                       
    eta3.jrule <- sqrt((exp(var3.jrule)-delta3.jrule)/delta3.jrule)                       
    gamma12.jrule <- eta1.jrule-eta2.jrule
    gamma13.jrule <- eta1.jrule-eta3.jrule
    gamma23.jrule <- eta2.jrule-eta3.jrule
    L.rBaye12_jrule[i] <-quantile(gamma12.jrule,alpha/2)                                  
    if(is.na(L.rBaye12_jrule[i])){next}
    if(is.nan(L.rBaye12_jrule[i])){next}
    U.rBaye12_jrule[i] <-quantile(gamma12.jrule,(1-alpha/2))                              
    if(is.na(U.rBaye12_jrule[i])){next}
    if(is.nan(U.rBaye12_jrule[i])){next}
    CP.rBaye12_jrule[i] <- ifelse(L.rBaye12_jrule[i]< gamma12 && gamma12 <U.rBaye12_jrule[i],1,0)  
    Length.rBaye12_jrule[i] <- U.rBaye12_jrule[i]-L.rBaye12_jrule[i]                               
    L.rBaye13_jrule[i] <-quantile(gamma13.jrule,alpha/2)                                  
    if(is.na(L.rBaye13_jrule[i])){next}
    if(is.nan(L.rBaye13_jrule[i])){next}
    U.rBaye13_jrule[i] <-quantile(gamma13.jrule,(1-alpha/2))                              
    if(is.na(U.rBaye13_jrule[i])){next}
    if(is.nan(U.rBaye13_jrule[i])){next}
    CP.rBaye13_jrule[i] <- ifelse(L.rBaye13_jrule[i]< gamma13 && gamma13 <U.rBaye13_jrule[i],1,0)  
    Length.rBaye13_jrule[i] <- U.rBaye13_jrule[i]-L.rBaye13_jrule[i]                               
    L.rBaye23_jrule[i] <-quantile(gamma23.jrule,alpha/2)                                  
    if(is.na(L.rBaye23_jrule[i])){next}
    if(is.nan(L.rBaye23_jrule[i])){next}
    U.rBaye23_jrule[i] <-quantile(gamma23.jrule,(1-alpha/2))                              
    if(is.na(U.rBaye23_jrule[i])){next}
    if(is.nan(U.rBaye23_jrule[i])){next}
    CP.rBaye23_jrule[i] <- ifelse(L.rBaye23_jrule[i]< gamma23 && gamma23 <U.rBaye23_jrule[i],1,0)     
    Length.rBaye23_jrule[i] <- U.rBaye23_jrule[i]-L.rBaye23_jrule[i]                                  
    CP.rBaye_jrule[i] <- (CP.rBaye12_jrule[i]+CP.rBaye13_jrule[i]+CP.rBaye23_jrule[i])/3              
    AL.rBaye_jrule[i] <- (Length.rBaye12_jrule[i]+Length.rBaye13_jrule[i]+Length.rBaye23_jrule[i])/3                                
    
    #Equitailed Bayesian based on uniform prior................................................................
    beta1.uni <- rbeta(m,size_nonzero1+1,(n1-size_nonzero1)+1)                        
    beta2.uni <- rbeta(m,size_nonzero2+1,(n2-size_nonzero2)+1)                        
    beta3.uni <- rbeta(m,size_nonzero3+1,(n3-size_nonzero3)+1)                        
    var1.uni <- rinvgamma(m,(size_nonzero1-2)/2,((size_nonzero1-2)*var_hat1)/2,scale = 1/rate)
    var2.uni <- rinvgamma(m,(size_nonzero2-2)/2,((size_nonzero2-2)*var_hat2)/2,scale = 1/rate)
    var3.uni <- rinvgamma(m,(size_nonzero3-2)/2,((size_nonzero3-2)*var_hat3)/2,scale = 1/rate)
    delta1.uni <- beta1.uni                                                           
    delta2.uni <- beta2.uni                                                           
    delta3.uni <- beta3.uni                                                           
    eta1.uni <-  sqrt((exp(var1.uni)-delta1.uni)/delta1.uni)                          
    eta2.uni <-  sqrt((exp(var2.uni)-delta2.uni)/delta2.uni)                          
    eta3.uni <-  sqrt((exp(var3.uni)-delta3.uni)/delta3.uni)                          
    gamma12.uni <- eta1.uni-eta2.uni
    gamma13.uni <- eta1.uni-eta3.uni
    gamma23.uni <- eta2.uni-eta3.uni
    L.rBaye12_uni[i] <-quantile(gamma12.uni,alpha/2)                                  
    if(is.na(L.rBaye12_uni[i])){next}
    if(is.nan(L.rBaye12_uni[i])){next}
    U.rBaye12_uni[i] <-quantile(gamma12.uni,(1-alpha/2))                              
    if(is.na(U.rBaye12_uni[i])){next}
    if(is.nan(U.rBaye12_uni[i])){next}
    CP.rBaye12_uni[i] <- ifelse(L.rBaye12_uni[i]< gamma12 && gamma12 <U.rBaye12_uni[i],1,0)    
    Length.rBaye12_uni[i] <- U.rBaye12_uni[i]-L.rBaye12_uni[i]                                 
    L.rBaye13_uni[i] <-quantile(gamma13.uni,alpha/2)                                  
    if(is.na(L.rBaye13_uni[i])){next}
    if(is.nan(L.rBaye13_uni[i])){next}
    U.rBaye13_uni[i] <-quantile(gamma13.uni,(1-alpha/2))                              
    if(is.na(U.rBaye13_uni[i])){next}
    if(is.nan(U.rBaye13_uni[i])){next}
    CP.rBaye13_uni[i] <- ifelse(L.rBaye13_uni[i]< gamma13 && gamma13 <U.rBaye13_uni[i],1,0)    
    Length.rBaye13_uni[i] <- U.rBaye13_uni[i]-L.rBaye13_uni[i]                                 
    L.rBaye23_uni[i] <-quantile(gamma23.uni,alpha/2)                                  
    if(is.na(L.rBaye23_uni[i])){next}
    if(is.nan(L.rBaye23_uni[i])){next}
    U.rBaye23_uni[i] <-quantile(gamma23.uni,(1-alpha/2))                              
    if(is.na(U.rBaye23_uni[i])){next}
    if(is.nan(U.rBaye23_uni[i])){next}
    CP.rBaye23_uni[i] <- ifelse(L.rBaye23_uni[i]< gamma23 && gamma23 <U.rBaye23_uni[i],1,0)    
    Length.rBaye23_uni[i] <- U.rBaye23_uni[i]-L.rBaye23_uni[i]                                 
    CP.rBaye_uni[i] <- (CP.rBaye12_uni[i]+CP.rBaye13_uni[i]+CP.rBaye23_uni[i])/3               
    AL.rBaye_uni[i] <- (Length.rBaye12_uni[i]+Length.rBaye13_uni[i]+Length.rBaye23_uni[i])/3                             
    
    #HPD Bayesian based on Jeffreys' rule prior..............................................................
    hpd12.jrule <- hdi(gamma12.jrule,0.95)                                            
    L.rBaye12_hpd.jrule[i] <- hpd12.jrule[1]                                          
    U.rBaye12_hpd.jrule[i] <- hpd12.jrule[2]                                          
    CP.rBaye12_hpd.jrule[i] <- ifelse(L.rBaye12_hpd.jrule[i]< gamma12 && gamma12 <U.rBaye12_hpd.jrule[i],1,0)    
    Length.rBaye12_hpd.jrule[i] <- U.rBaye12_hpd.jrule[i]-L.rBaye12_hpd.jrule[i]                                 
    hpd13.jrule <- hdi(gamma13.jrule,0.95)                                            
    L.rBaye13_hpd.jrule[i] <- hpd13.jrule[1]                                          
    U.rBaye13_hpd.jrule[i] <- hpd13.jrule[2]                                          
    CP.rBaye13_hpd.jrule[i] <- ifelse(L.rBaye13_hpd.jrule[i]< gamma13 && gamma13 <U.rBaye13_hpd.jrule[i],1,0)    
    Length.rBaye13_hpd.jrule[i] <- U.rBaye13_hpd.jrule[i]-L.rBaye13_hpd.jrule[i]                                 
    hpd23.jrule <- hdi(gamma23.jrule,0.95)                                            
    L.rBaye23_hpd.jrule[i] <- hpd23.jrule[1]                                          
    U.rBaye23_hpd.jrule[i] <- hpd23.jrule[2]                                          
    CP.rBaye23_hpd.jrule[i] <- ifelse(L.rBaye23_hpd.jrule[i]< gamma23 && gamma23 <U.rBaye23_hpd.jrule[i],1,0)    
    Length.rBaye23_hpd.jrule[i] <- U.rBaye23_hpd.jrule[i]-L.rBaye23_hpd.jrule[i]                                 
    CP.rBaye_hpd.jrule[i] <- (CP.rBaye12_hpd.jrule[i]+CP.rBaye13_hpd.jrule[i]+CP.rBaye23_hpd.jrule[i])/3               
    AL.rBaye_hpd.jrule[i] <- (Length.rBaye12_hpd.jrule[i]+Length.rBaye13_hpd.jrule[i]+Length.rBaye23_hpd.jrule[i])/3
    
    #HPD Bayesian based on Uniform prior.......................................................................
    hpd12.uni <- hdi(gamma12.uni,0.95)                                                
    L.rBaye12_hpd.uni[i] <- hpd12.uni[1]                                              
    U.rBaye12_hpd.uni[i] <- hpd12.uni[2]                                              
    CP.rBaye12_hpd.uni[i] <- ifelse(L.rBaye12_hpd.uni[i]< gamma12 && gamma12 <U.rBaye12_hpd.uni[i],1,0)    
    Length.rBaye12_hpd.uni[i] <- U.rBaye12_hpd.uni[i]-L.rBaye12_hpd.uni[i]                                 
    hpd13.uni <- hdi(gamma13.uni,0.95)                                                
    L.rBaye13_hpd.uni[i] <- hpd13.uni[1]                                              
    U.rBaye13_hpd.uni[i] <- hpd13.uni[2]                                              
    CP.rBaye13_hpd.uni[i] <- ifelse(L.rBaye13_hpd.uni[i]< gamma13 && gamma13 <U.rBaye13_hpd.uni[i],1,0)    
    Length.rBaye13_hpd.uni[i] <- U.rBaye13_hpd.uni[i]-L.rBaye13_hpd.uni[i]                                 
    hpd23.uni <- hdi(gamma23.uni,0.95)                                                
    L.rBaye23_hpd.uni[i] <- hpd23.uni[1]                                              
    U.rBaye23_hpd.uni[i] <- hpd23.uni[2]                                              
    CP.rBaye23_hpd.uni[i] <- ifelse(L.rBaye23_hpd.uni[i]< gamma23 && gamma23 <U.rBaye23_hpd.uni[i],1,0)        
    Length.rBaye23_hpd.uni[i] <- U.rBaye23_hpd.uni[i]-L.rBaye23_hpd.uni[i]                                     
    CP.rBaye_hpd.uni[i] <- (CP.rBaye12_hpd.uni[i]+CP.rBaye13_hpd.uni[i]+CP.rBaye23_hpd.uni[i])/3               
    AL.rBaye_hpd.uni[i] <- (Length.rBaye12_hpd.uni[i]+Length.rBaye13_hpd.uni[i]+Length.rBaye23_hpd.uni[i])/3
    
    #MOVER based on Wilson score method........................................................................
    s <- qnorm(1-(alpha/2),0,1)
    
    df1 <- size_nonzero1-1
    bound_U1 <- qchisq(alpha/2,df1) 
    bound_L1 <- qchisq(1-(alpha/2),df1) 
    L_Var1 <- ((size_nonzero1-1)*var_hat1)/bound_L1                       
    U_Var1 <- ((size_nonzero1-1)*var_hat1)/bound_U1
    
    df2 <- size_nonzero2-1
    bound_U2 <- qchisq(alpha/2,df2) 
    bound_L2 <- qchisq(1-(alpha/2),df2) 
    L_Var2 <- ((size_nonzero2-1)*var_hat2)/bound_L2                       
    U_Var2 <- ((size_nonzero2-1)*var_hat2)/bound_U2
    
    df3 <- size_nonzero3-1
    bound_U3 <- qchisq(alpha/2,df3) 
    bound_L3 <- qchisq(1-(alpha/2),df3) 
    L_Var3 <- ((size_nonzero3-1)*var_hat3)/bound_L3                       
    U_Var3 <- ((size_nonzero3-1)*var_hat3)/bound_U3
    
    deltahat_W1 <- (size_nonzero1+((s^2)/2))/(n1+(s^2))
    var_deltahat_W1 <-(((n1-size_nonzero1)*size_nonzero1/n1)+((s^2)/4))
    L_delta1_W <- deltahat_W1-(s*sqrt(var_deltahat_W1)/(n1+(s^2)))
    U_delta1_W <- deltahat_W1+(s*sqrt(var_deltahat_W1)/(n1+(s^2)))
    
    deltahat_W2 <- (size_nonzero2+((s^2)/2))/(n2+(s^2))
    var_deltahat_W2 <-(((n2-size_nonzero2)*size_nonzero2/n2)+((s^2)/4))
    L_delta2_W <- deltahat_W2-(s*sqrt(var_deltahat_W2)/(n2+(s^2)))
    U_delta2_W <- deltahat_W2+(s*sqrt(var_deltahat_W2)/(n2+(s^2)))
    
    deltahat_W3 <- (size_nonzero3+((s^2)/2))/(n3+(s^2))
    var_deltahat_W3 <-(((n3-size_nonzero3)*size_nonzero3/n3)+((s^2)/4))
    L_delta3_W <- deltahat_W3-(s*sqrt(var_deltahat_W3)/(n3+(s^2)))
    U_delta3_W <- deltahat_W3+(s*sqrt(var_deltahat_W3)/(n3+(s^2)))
    
    l1_W <- sqrt((exp(L_Var1)-U_delta1_W)/U_delta1_W)
    u1_W <- sqrt((exp(U_Var1)-L_delta1_W)/L_delta1_W)
    l2_W <- sqrt((exp(L_Var2)-U_delta2_W)/U_delta2_W)
    u2_W <- sqrt((exp(U_Var2)-L_delta2_W)/L_delta2_W)
    l3_W <- sqrt((exp(L_Var3)-U_delta3_W)/U_delta3_W)
    u3_W <- sqrt((exp(U_Var3)-L_delta3_W)/L_delta3_W) 
    
    L.rMOVER12_W[i] <-(eta_hat1-eta_hat2)-sqrt(((eta_hat1-l1_W)^2)+((u2_W-eta_hat2)^2))           
    U.rMOVER12_W[i] <-(eta_hat1-eta_hat2)+sqrt(((u1_W-eta_hat1)^2)+((eta_hat2-l2_W)^2))           
    CP.rMOVER12_W[i] <- ifelse(L.rMOVER12_W[i]< gamma12 && gamma12 <U.rMOVER12_W[i],1,0)          
    Length.rMOVER12_W[i] <- U.rMOVER12_W[i]-L.rMOVER12_W[i]                                       
    L.rMOVER13_W[i] <-(eta_hat1-eta_hat3)-sqrt(((eta_hat1-l1_W)^2)+((u3_W-eta_hat3)^2))           
    U.rMOVER13_W[i] <-(eta_hat1-eta_hat3)+sqrt(((u1_W-eta_hat1)^2)+((eta_hat3-l3_W)^2))           
    CP.rMOVER13_W[i] <- ifelse(L.rMOVER13_W[i]< gamma13 && gamma13 <U.rMOVER13_W[i],1,0)          
    Length.rMOVER13_W[i] <- U.rMOVER13_W[i]-L.rMOVER13_W[i]                                       
    L.rMOVER23_W[i] <-(eta_hat2-eta_hat3)-sqrt(((eta_hat2-l2_W)^2)+((u3_W-eta_hat3)^2))           
    U.rMOVER23_W[i] <-(eta_hat2-eta_hat3)+sqrt(((u2_W-eta_hat2)^2)+((eta_hat3-l3_W)^2))           
    CP.rMOVER23_W[i] <- ifelse(L.rMOVER23_W[i]< gamma23 && gamma23 <U.rMOVER23_W[i],1,0)          
    Length.rMOVER23_W[i] <- U.rMOVER23_W[i]-L.rMOVER23_W[i]                                       
    CP.rMOVER_W[i] <- (CP.rMOVER12_W[i]+CP.rMOVER13_W[i]+CP.rMOVER23_W[i])/3                      
    AL.rMOVER_W[i] <- (Length.rMOVER12_W[i]+Length.rMOVER13_W[i]+Length.rMOVER23_W[i])/3  
    
    i=i+1
    print(i)
  }                                                                         #end loop i
  
  ACP.rFGCI <- mean(CP.rFGCI)                              
  Ex.length.rFGCI <- mean(AL.rFGCI)                      
  ACP.rBaye_jrule <- mean(CP.rBaye_jrule)                  
  Ex.length.rBaye_jrule <- mean(AL.rBaye_jrule)          
  ACP.rBaye_uni <- mean(CP.rBaye_uni)                      
  Ex.length.rBaye_uni <- mean(AL.rBaye_uni)               
  ACP.rBaye_hpd.jrule <- mean(CP.rBaye_hpd.jrule)          
  Ex.length.rBaye_hpd.jrule <- mean(AL.rBaye_hpd.jrule)   
  ACP.rBaye_hpd.uni <- mean(CP.rBaye_hpd.uni)              
  Ex.length.rBaye_hpd.uni <- mean(AL.rBaye_hpd.uni)      
  ACP.rMOVER_W <- mean(CP.rMOVER_W)                        
  Ex.length.rMOVER_W <- mean(AL.rMOVER_W)                 
  elapseTime = proc.time()-starttime
  cat("M =",M,"m =",m,"n1 =",n1,"n2 =",n2,"n3 =",n3,"delta1 =",delta1,"delta2 =",delta2,
          "delta3 =",delta3,"var1 =",var1,"var2 =",var2,"var3 =",var3,"\n",
      "E1(x_0)=",n1*delta1,"E2(x_0)=",n2*delta2,"E3(x_0)=",n3*delta3,"cv1=",cv1,"cv2=",cv2,"cv3=",cv3,
          "gamma12",gamma12,"gamma13",gamma13,"gamma23",gamma23,"\n",
      "CP|rFGCI =",ACP.rFGCI,"Length|rFGCI =",Ex.length.rFGCI,"\n",
      "CP|rBaye_jrule =",ACP.rBaye_jrule,"Length|rBaye_jrule =",Ex.length.rBaye_jrule,"\n",
      "CP|rBaye_uni =",ACP.rBaye_uni,"Length|rBaye_uni =",Ex.length.rBaye_uni,"\n",
      "CP|rBaye_hpd.jrule =",ACP.rBaye_hpd.jrule,"Length|rBaye_hpd.jrule =",Ex.length.rBaye_hpd.jrule,"\n",
      "CP|rBaye_hpd.uni =",ACP.rBaye_hpd.uni,"Length|rBaye_hpd.uni =",Ex.length.rBaye_hpd.uni,"\n",
      "CP|rMOVER_W =",ACP.rMOVER_W,"Length|rMOVER_W =",Ex.length.rMOVER_W)
  result <- c("M =",M,"m =",m,"n1 =",n1,"n2 =",n2,"n3 =",n3,"delta1 =",delta1,"delta2 =",delta2,
                  "delta3 =",delta3,"var1 =",var1,"var2 =",var2,"var3 =",var3,"\n",
              "E1(x_0)=",n1*delta1,"E2(x_0)=",n2*delta2,"E3(x_0)=",n3*delta3,"cv1=",cv1,"cv2=",cv2,"cv3=",cv3,
                  "gamma12",gamma12,"gamma13",gamma13,"gamma23",gamma23,"\n",
              "CP|rFGCI =",ACP.rFGCI,"Length|rFGCI =",Ex.length.rFGCI,"\n",
              "CP|rBaye_jrule =",ACP.rBaye_jrule,"Length|rBaye_jrule =",Ex.length.rBaye_jrule,"\n",
              "CP|rBaye_uni =",ACP.rBaye_uni,"Length|rBaye_uni =",Ex.length.rBaye_uni,"\n",
              "CP|rBaye_hpd.jrule =",ACP.rBaye_hpd.jrule,"Length|rBaye_hpd.jrule =",Ex.length.rBaye_hpd.jrule,"\n",
              "CP|rBaye_hpd.uni =",ACP.rBaye_hpd.uni,"Length|rBaye_hpd.uni =",Ex.length.rBaye_hpd.uni,"\n",
              "CP|rMOVER_W =",ACP.rMOVER_W,"Length|rMOVER_W =",Ex.length.rMOVER_W,"\n",
              "Time =",elapseTime["elapsed"])
  return(result)  
}