##################################################################################
# R code to construct SCIs for all possible ratios of variances of several ZILN 
# based on Bayesian GPQ and PB approaches with        
# Daily precipitation Datasets 
########################################S##########################################


    ##################################################################################
    ##################################################################################
    # Three sample groups
    # k = 3
    ##################################################################################
    ##################################################################################

    #Simulation runs
    M = 5000
    m = 2500
    
    
    # Function for generating random sample drawn from delta-lognormal distributions with k-th populations
    library(EnvStats)
    
    xGen = function(k,nTotal,Ex_LN,cv_LN,deltaTotal){
      sum_nonzero_value2 = rep(0,k)
      mu_hat = rep(0,k)
      sigmaSq_hat= rep(0,k)
      delta_hatnon= rep(0,k)
      size_non = rep(0,k)
      size_zero = rep(0,k)
      xlist = list()
      
      
      for (l in 1:k){
        #print(l)
        while(TRUE){
          x = rzmlnormAlt(nTotal[l], Ex_LN[l], cv_LN[l], deltaTotal[l])
          xlist[l] = list(x)
          
          Temp.nonzero_value = log(x[which(x!=0)])                         
          Temp.size_nonzero =length(Temp.nonzero_value)                          
          if(Temp.size_nonzero < (nTotal[l]-nTotal[l]*deltaTotal[l])){
            next
          }
          nonzero_value = Temp.nonzero_value
          size_nonzero = Temp.size_nonzero
          sum_nonzero_value2[l] = sum(nonzero_value^2)
          mu_hat[l] = mean(nonzero_value)
          sigmaSq_hat[l] = var(nonzero_value)
          
          delta_hatnon[l] = size_nonzero/nTotal[l]
          size_non[l] = size_nonzero
          size_zero[l] = nTotal[l]-size_nonzero
          
          break
        }
      }
      return(list("data" = cbind(size_non, size_zero, mu_hat, sigmaSq_hat, delta_hatnon, sum_nonzero_value2),"xlist"=xlist))
    }
    
    
    
    # Parameter estimations
    writeDataToFile=function(x,k,rep,currentFolderPath){
      datatable = t(data.frame(x))
      filename = paste(currentFolderPath,'/D',k,'_xlist.csv',sep="")
      write.table( datatable, filename  , append= T, sep=',',col.names=FALSE,row.names = FALSE )
    }
    
    
    createDatasource = function(nCase,sCase,dCase){
      
      ###############################Delta Varied########################################
      folderSourceName = paste("n",nCase,"s",sCase,"d",dCase, sep = "")
      print(folderSourceName)
      ###################################################################################
      #Create folder
      dir.create(file.path("./DataSource/"), showWarnings = FALSE)
      
      ## set current folder to keep all data  ./Datasource/Data_m_2500_n1_100n2_100.......
      currentFolderPath = paste("./DataSource/",folderSourceName,sep="")
      dir.create(file.path(currentFolderPath),showWarnings = FALSE)
      
      return (currentFolderPath)
      
    }
    
      ntotal_allcase = matrix(c( 30,   30,    30,
                                 50,   50,    50,
                                 100,  100,   100,
                                 30,   50,    100,  
                                 50,   100,   200,
                                 100,  100,   200), ncol = 3,byrow = TRUE)
      
      sigmaSq_allcase = matrix(c(3,  5,    7), ncol = 3,byrow = TRUE) 
      
      delta_allcase = matrix(c( 0.1,  0.2,  0.3,
                                0.1,  0.3,  0.5,
                                0.3,  0.5,  0.5), ncol = 3,byrow = TRUE) 
    
    
    for (id_ntotal_alcase in seq(1:nrow(ntotal_allcase))){
      config_ntotal =ntotal_allcase[id_ntotal_alcase,]
      
      for (id_sigmaSq_allcase in seq(1:nrow(sigmaSq_allcase))){
        config_sigmaSq =sigmaSq_allcase[id_sigmaSq_allcase,]
        
        for (id_delta_allcase in seq(1:nrow(delta_allcase))){
          config_delta =delta_allcase[id_delta_allcase,]
          
          currentFolderPath = createDatasource(id_ntotal_alcase,id_sigmaSq_allcase,id_delta_allcase)
          
          allData = data.frame(data.frame(config_ntotal),data.frame(config_sigmaSq),data.frame(config_delta))
          write.csv(allData,file=paste(currentFolderPath,"/config.csv",sep=""),row.names = FALSE)
        }
        
      }
      
    }
    
    listFolder = list.dirs("./DataSource/", recursive=FALSE)
    
    eachFolder = listFolder[1]
    print(listFolder)
    for (eachFolder in listFolder){
      currentFolder = eachFolder
      configfile = read.csv(file=paste(currentFolder,"/config.csv",sep = ""))
      nTotal = configfile$config_ntotal
      sigmaSqTotal = configfile$config_sigmaSq
      deltaTotal = configfile$config_delta
      
      
      k = length(nTotal)
      
      muTotal =  -sigmaSqTotal/2
      mu1 = muTotal[1]
      mu2 = muTotal[2]
      mu3 = muTotal[3]
      
      deltanonTotal = 1 - deltaTotal
      size_nonzeroTotal = nTotal*(1-deltaTotal)
      
      
      Ex_LN = exp(muTotal+(sigmaSqTotal/2))
      cv_LN = sqrt(exp(sigmaSqTotal)-1)
      theta = log(deltanonTotal)+(2*muTotal+sigmaSqTotal)+log(exp(sigmaSqTotal)-deltanonTotal)
      theta12 = theta[1]-theta[2]
      theta13 = theta[1]-theta[3]
      theta23 = theta[2]-theta[3]
      
      
      
      AllDataset = data.frame()
      ## Check file .xlxs is exist it will be deleted
      
      delfiles = dir(path=currentFolder ,pattern="D.*xlist.csv")
      # Check its existence
      file.remove(file.path(currentFolder, delfiles))
      
      idx_rep =1
      for ( idx_rep in seq(1:M)){
        #Start generate data
        packData = xGen(k, nTotal, Ex_LN, cv_LN, deltaTotal)
        dataset= data.frame(loop = rep(idx_rep,k),num_k= 1:k)
        data  = packData$data
        xlist = packData$xlist
        size_nonzero = data
        size_zero = data[,2]
        mu_hat = data[,3]
        sigmaSq_hat  = data[,4]
        delta_hatnon = data[,5]
        delta_hat  = 1-delta_hatnon
        sum_nonzero_value2 = data[,6]
        dataset = cbind(dataset, data)
        
        ## Store Dataset each loop
        AllDataset = rbind(AllDataset,dataset)
        
        
        # Write raw xlist data
        for (idx_xlist in seq(1:length(xlist))){
          eachlist = xlist[idx_xlist]
          writeDataToFile(eachlist,idx_xlist,idx_rep,currentFolder)
        }
        
      }
      write.table( AllDataset, paste(currentFolder,"/basicSummary_K",k,".csv",sep="")  , sep=',',row.names = FALSE )
    }
    



    ##############################################################################
    #                 Algorithm 1: B-PMB                                         #
    ##############################################################################   
    
    library(rstan)            #Stan file
    
    ### This function for loading the raw dataset
    loadXlist = function(folderselect,Kselect){
      D1 = read.csv(file=paste(folderselect,"/D",Kselect,"_xlist.csv",sep = ""),header = FALSE)
      return (D1)
    }
    
    
    resultoutput = "./BPM_result/"
    
    #Delete all result
    delfiles = dir(path=resultoutput ,pattern="*.txt")
    file.remove(file.path(resultoutput, delfiles))
    
    M = 5000
    m = 2500
    
    listFolder = list.dirs("./DataSource/", recursive=FALSE)
    
    eachFolder = listFolder[1]
    print(listFolder)
    for (eachFolder in listFolder){
      
      print(paste("Start Folder",eachFolder,sep=""))
      ## read data from config file
      currentFolder = eachFolder
      configfile = read.csv(file=paste(currentFolder,"/config.csv",sep = ""))
      nTotal = configfile$config_ntotal
      sigmaSqTotal = configfile$config_sigmaSq
      deltaTotal = configfile$config_delta
      muTotal = -sigmaSqTotal/2
      k = length(nTotal)
      
      
      # Calculate parameter from data in config file
      deltanonTotal = 1 - deltaTotal
      size_nonzeroTotal = nTotal*(1-deltaTotal)
      Ex_LN = exp(muTotal+(sigmaSqTotal/2))
      cv_LN = sqrt(exp(sigmaSqTotal)-1)
      theta = log(deltanonTotal)+(2*muTotal+sigmaSqTotal)+log(exp(sigmaSqTotal)-deltanonTotal)
      theta12 = theta[1]-theta[2]
      theta13 = theta[1]-theta[3]
      theta23 = theta[2]-theta[3]
      
      # Read summary data
      summaryData = read.csv(file= paste(currentFolder,"/basicSummary_K",k,".csv",sep = ""))
      
      # B-PMB: Bayesian-based probability-matching-beta prior 
      CI.BPM = data.frame()
      D.BPM  = rep(0,m)
      
      rep = 1
      dataframe.x1 = loadXlist(currentFolder,1)
      dataframe.x2 = loadXlist(currentFolder,2)
      dataframe.x3 = loadXlist(currentFolder,3)
      
      for(rep in seq(1:M)){
        
        #Select data from repeat loop number
        currentData = summaryData[which(summaryData$loop==rep),]
        
        #Extract each parameter 
        delta_hatnon = currentData$delta_hatnon
        delta_hat = 1-delta_hatnon
        size_nonzero = currentData$size_non
        size_zero = currentData$size_zero
        mu_hat = currentData$mu_hat
        sigmaSq_hat = currentData$sigmaSq_hat
        n = size_nonzero + size_zero
        
        thetahat = log(delta_hatnon) + (2*(mu_hat+sigmaSq_hat))
        thetahat.12 = thetahat[1]-thetahat[2]
        thetahat.13 = thetahat[1]-thetahat[3]
        thetahat.23 = thetahat[2]-thetahat[3]
        
        Var_varDLN1 = delta_hat/(n*delta_hatnon)
        Var_varDLN2 = (sigmaSq_hat/size_nonzero) + (2*sigmaSq_hat^2/(size_nonzero-1))
        Var_varDLN = Var_varDLN1 + (4*Var_varDLN2)
        
        Var_varDLN12 = Var_varDLN[1]+Var_varDLN[2]
        Var_varDLN13 = Var_varDLN[1]+Var_varDLN[3]
        Var_varDLN23 = Var_varDLN[2]+Var_varDLN[3]
        
        
       ############################### BPM ######################################## 
        ## chain = 4, iter = 1250 to obtain m = 4*(1250/2) = 2500  
        ## refresh = 0 denoted as Stan running without output details
        
        x1 = as.numeric(dataframe.x1[rep,])
        x2 = as.numeric(dataframe.x2[rep,])
        x3 = as.numeric(dataframe.x3[rep,])
        
        
        fit.pm = stan('Stan_PMBprior_K3.stan', data=list(n1=n[1], sizenonzero1=size_nonzero[1], 
                                                                 y1=x1, sizezero1=size_zero[1], mu1=mu_hat[1], 
                                                                 sigmaSq1=sigmaSq_hat[1], delta1=delta_hat[1],
                                                                 
                                                                 n2=n[2], sizenonzero2=size_nonzero[2], 
                                                                 y2=x2, sizezero2=size_zero[2], mu2=mu_hat[2], 
                                                                 sigmaSq2=sigmaSq_hat[2], delta2=delta_hat[2],
                                                                 
                                                                 n3=n[3], sizenonzero3=size_nonzero[3], 
                                                                 y3=x3, sizezero3=size_zero[3], mu3=mu_hat[3], 
                                                                 sigmaSq3=sigmaSq_hat[3], delta3=delta_hat[3]), 
                      iter = 1250, refresh = 0)
        
        # Extracting the posterior draws
        post.mu.pm1 = extract(fit.pm)$mu1
        post.sigmaSq.pm1 = extract(fit.pm)$sigmaSq1
        post.delta.pm1 = extract(fit.pm)$delta1
        post_varDLN.pm1 = log(1-post.delta.pm1)+2*(post.mu.pm1+post.sigmaSq.pm1)
        
        post.mu.pm2 = extract(fit.pm)$mu2
        post.sigmaSq.pm2 = extract(fit.pm)$sigmaSq2
        post.delta.pm2 = extract(fit.pm)$delta2
        post_varDLN.pm2 = log(1-post.delta.pm2)+2*(post.mu.pm2+post.sigmaSq.pm2)
        
        post.mu.pm3 = extract(fit.pm)$mu3
        post.sigmaSq.pm3 = extract(fit.pm)$sigmaSq3
        post.delta.pm3 = extract(fit.pm)$delta3
        post_varDLN.pm3 = log(1-post.delta.pm3)+2*(post.mu.pm3+post.sigmaSq.pm3)
        
        post_varDLN.pm12 = post_varDLN.pm1-post_varDLN.pm2
        post_varDLN.pm13 = post_varDLN.pm1-post_varDLN.pm3
        post_varDLN.pm23 = post_varDLN.pm2-post_varDLN.pm3
        
        for(idx.BPM in 1:m){
          D12.BPM = abs(thetahat.12 - post_varDLN.pm12[idx.BPM])
          D13.BPM = abs(thetahat.13 - post_varDLN.pm13[idx.BPM])
          D23.BPM = abs(thetahat.23 - post_varDLN.pm23[idx.BPM])
          
          A.BPM = max(D12.BPM, D13.BPM, D23.BPM)
          B.BPM = min(D12.BPM, D13.BPM, D23.BPM)
          
          D.BPM[idx.BPM] = A.BPM - B.BPM
          
        }
        
        
        q.BPM =   quantile(D.BPM,0.95,type=8)
        L.BPM12 = thetahat.12 - q.BPM                     
        U.BPM12 = thetahat.12 + q.BPM
        CI.BPM12 = data.frame(rep,"CI.BPM12",L.BPM12,U.BPM12,theta12)
        names(CI.BPM12) = c("rep","CI","lower","upper","theta")
        CI.BPM = rbind(CI.BPM,CI.BPM12) 
        
        L.BPM13 = thetahat.13 - q.BPM                       
        U.BPM13 = thetahat.13 + q.BPM
        CI.BPM13 = data.frame(rep,"CI.BPM13",L.BPM13,U.BPM13,theta13)
        names(CI.BPM13) = c("rep","CI","lower","upper","theta")
        CI.BPM = rbind(CI.BPM,CI.BPM13) 
        
        L.BPM23 = thetahat.23 - q.BPM                    
        U.BPM23 = thetahat.23 + q.BPM
        CI.BPM23 = data.frame(rep,"CI.BPM23",L.BPM23,U.BPM23,theta23)
        names(CI.BPM23) = c("rep","CI","lower","upper","theta")
        CI.BPM = rbind(CI.BPM,CI.BPM23) 
        
        
        if(rep%%100 ==0){
          print(rep)
        }
      }
      
      CI.BPM['CP.LE'] = 0
      CI.BPM['CP'] = 0
      CI.BPM['CP.UE'] = 0
      CI.BPM[which(CI.BPM["lower"]<CI.BPM["theta"] & CI.BPM["upper"]>CI.BPM["theta"]),"CP"]=1
      CI.BPM[which(CI.BPM["lower"]>CI.BPM["theta"]),"CP.LE"]=1
      CI.BPM[which(CI.BPM["upper"]<CI.BPM["theta"]),"CP.UE"]=1
      CI.BPM["Length"] = CI.BPM["upper"] - CI.BPM["lower"]
      ACI.BPM = colMeans(CI.BPM[,c("CP.LE","CP","CP.UE","Length")])
      ACP.BPM = as.numeric(ACI.BPM["CP"])*100
      ACP.LE.BPM = as.numeric(ACI.BPM["CP.LE"])*100
      ACP.UE.BPM = as.numeric(ACI.BPM["CP.UE"])*100
      Ex.length.BPM = as.numeric(ACI.BPM["Length"])	
      
      str1 = paste("M=",M,"m=",m, "mu=",paste(muTotal,collapse = ' '),"k=",k, "nTotal=",paste(nTotal,collapse = ' '))
      str2 = paste("sigmaSqTotal=", paste(sigmaSqTotal,collapse = ' '),"deltaTotal=",paste(deltaTotal,collapse = ' '),"theta=",paste(theta,collapse = ' '))   
      
      str3 =paste( "LE|BPM =",paste(ACP.LE.BPM,collapse = ' '),"CP|BPM =", paste(ACP.BPM,collapse = ' '),"UE|BPM =", paste(ACP.UE.BPM,collapse = ' '), "AL|BPM =", paste(Ex.length.BPM,collapse = ' '))
      
      alloutput = paste(str1,str2,str3,collapse = ' ')
      folderName = unlist(strsplit(currentFolder, "//"))
      dir.create(file.path(resultoutput), showWarnings = FALSE)
      write(alloutput, file = paste(resultoutput,"/Results_BPM.txt",sep = ""), append = TRUE)
    }


    ##############################################################################
    #                 Algorithm 2 : B-RB                                        #
    ##############################################################################   
    library(rstan)            #To use for read stan file
    
    ### This function for loading the raw dataset
    loadXlist = function(folderselect,Kselect){
      D1 = read.csv(file=paste(folderselect,"/D",Kselect,"_xlist.csv",sep = ""),header = FALSE)
      return (D1)
    }
    
    resultoutput = "./BR_result/"
    
    #Delete all result
    delfiles = dir(path=resultoutput ,pattern="*.txt")
    file.remove(file.path(resultoutput, delfiles))
    
    M = 5000
    m = 2500
    
    
    listFolder = list.dirs("./DataSource/", recursive=FALSE)
    
    eachFolder = listFolder[1]
    print(listFolder)
    for (eachFolder in listFolder){
      
      print(paste("Start Folder",eachFolder,sep=""))
      ## read data from config file
      currentFolder = eachFolder
      configfile = read.csv(file=paste(currentFolder,"/config.csv",sep = ""))
      nTotal = configfile$config_ntotal
      sigmaSqTotal = configfile$config_sigmaSq
      deltaTotal = configfile$config_delta
      muTotal = -sigmaSqTotal/2
      k = length(nTotal)
      
      
      # calculate parameter from data in config file
      deltanonTotal = 1 - deltaTotal
      size_nonzeroTotal = nTotal*(1-deltaTotal)
      Ex_LN = exp(muTotal+(sigmaSqTotal/2))
      cv_LN = sqrt(exp(sigmaSqTotal)-1)
      theta = log(deltanonTotal)+(2*muTotal+sigmaSqTotal)+log(exp(sigmaSqTotal)-deltanonTotal)
      theta12 = theta[1]-theta[2]
      theta13 = theta[1]-theta[3]
      theta23 = theta[2]-theta[3]
      
      #read summary data
      summaryData = read.csv(file= paste(currentFolder,"/basicSummary_K",k,".csv",sep = ""))
      
      #BR: Bayesian-based reference prior 
      CI.BR = data.frame()
      D.BR  = rep(0,m)
      
      rep = 1
      dataframe.x1 = loadXlist(currentFolder,1)
      dataframe.x2 = loadXlist(currentFolder,2)
      dataframe.x3 = loadXlist(currentFolder,3)
      for(rep in seq(1:M)){
        
        #Select data from repeat loop number
        currentData = summaryData[which(summaryData$loop==rep),]
        
        #Extract each parameter 
        delta_hatnon = currentData$delta_hatnon
        delta_hat = 1-delta_hatnon
        size_nonzero = currentData$size_non
        size_zero = currentData$size_zero
        mu_hat = currentData$mu_hat
        sigmaSq_hat = currentData$sigmaSq_hat
        n = size_nonzero + size_zero
        
        thetahat = log(delta_hatnon) + (2*(mu_hat+sigmaSq_hat))
        thetahat.12 = thetahat[1]-thetahat[2]
        thetahat.13 = thetahat[1]-thetahat[3]
        thetahat.23 = thetahat[2]-thetahat[3]
        
        Var_varDLN1 = delta_hat/(n*delta_hatnon)
        Var_varDLN2 = (sigmaSq_hat/size_nonzero) + (2*sigmaSq_hat^2/(size_nonzero-1))
        Var_varDLN = Var_varDLN1 + (4*Var_varDLN2)
        
        Var_varDLN12 = Var_varDLN[1]+Var_varDLN[2]
        Var_varDLN13 = Var_varDLN[1]+Var_varDLN[3]
        Var_varDLN23 = Var_varDLN[2]+Var_varDLN[3]
        
        
        
        ############################### BR ######################################## 
        ## chain = 4, iter = 1250 to obtain m = 4*(1250/2) = 2500  
        ## refresh = 0 denoted as Stan running without output details
        x1 = as.numeric(dataframe.x1[rep,])
        x2 = as.numeric(dataframe.x2[rep,])
        x3 = as.numeric(dataframe.x3[rep,])
        
        fit.ref = stan('Stan_RBprior_K3.stan', data=list(n1=n[1], sizenonzero1=size_nonzero[1], 
                                                               y1=x1, sizezero1=size_zero[1], mu1=mu_hat[1], 
                                                               sigmaSq1=sigmaSq_hat[1], delta1=delta_hat[1],
                                                               
                                                               n2=n[2], sizenonzero2=size_nonzero[2], 
                                                               y2=x2, sizezero2=size_zero[2], mu2=mu_hat[2], 
                                                               sigmaSq2=sigmaSq_hat[2], delta2=delta_hat[2],
                                                               
                                                               n3=n[3], sizenonzero3=size_nonzero[3], 
                                                               y3=x3, sizezero3=size_zero[3], mu3=mu_hat[3], 
                                                               sigmaSq3=sigmaSq_hat[3], delta3=delta_hat[3]),
                       iter = 1250, refresh = 0)
        
        ## Extracting the posterior draws
        post.mu.ref1 = extract(fit.ref)$mu1
        post.sigmaSq.ref1 = extract(fit.ref)$sigmaSq1
        post.delta.ref1 = extract(fit.ref)$delta1
        post_varDLN.ref1 = log(1-post.delta.ref1)+2*(post.mu.ref1+post.sigmaSq.ref1)
        
        post.mu.ref2 = extract(fit.ref)$mu2
        post.sigmaSq.ref2 = extract(fit.ref)$sigmaSq2
        post.delta.ref2 = extract(fit.ref)$delta2
        post_varDLN.ref2 = log(1-post.delta.ref2)+2*(post.mu.ref2+post.sigmaSq.ref2)
        
        post.mu.ref3 = extract(fit.ref)$mu3
        post.sigmaSq.ref3 = extract(fit.ref)$sigmaSq3
        post.delta.ref3 = extract(fit.ref)$delta3
        post_varDLN.ref3 = log(1-post.delta.ref3)+2*(post.mu.ref3+post.sigmaSq.ref3)
        
        post_varDLN.ref12 = post_varDLN.ref1-post_varDLN.ref2
        post_varDLN.ref13 = post_varDLN.ref1-post_varDLN.ref3
        post_varDLN.ref23 = post_varDLN.ref2-post_varDLN.ref3
        
        for(idx.BR in 1:m){
          D12.BR = abs(thetahat.12 - post_varDLN.ref12[idx.BR])
          D13.BR = abs(thetahat.13 - post_varDLN.ref13[idx.BR])
          D23.BR = abs(thetahat.23 - post_varDLN.ref23[idx.BR])
          
          A.BR = max(D12.BR, D13.BR, D23.BR)
          B.BR = min(D12.BR, D13.BR, D23.BR)
          
          D.BR[idx.BR] = A.BR - B.BR
          
        }
        
        
        q.BR =   quantile(D.BR,0.95,type=8)
        L.BR12 = thetahat.12 - q.BR                     
        U.BR12 = thetahat.12 + q.BR
        CI.BR12 = data.frame(rep,"CI.BR12",L.BR12,U.BR12,theta12)
        names(CI.BR12) = c("rep","CI","lower","upper","theta")
        CI.BR = rbind(CI.BR,CI.BR12) 
        
        L.BR13 = thetahat.13 - q.BR                       
        U.BR13 = thetahat.13 + q.BR
        CI.BR13 = data.frame(rep,"CI.BR13",L.BR13,U.BR13,theta13)
        names(CI.BR13) = c("rep","CI","lower","upper","theta")
        CI.BR = rbind(CI.BR,CI.BR13) 
        
        L.BR23 = thetahat.23 - q.BR                    
        U.BR23 = thetahat.23 + q.BR
        CI.BR23 = data.frame(rep,"CI.BR23",L.BR23,U.BR23,theta23)
        names(CI.BR23) = c("rep","CI","lower","upper","theta")
        CI.BR = rbind(CI.BR,CI.BR23) 
        
        
        if(rep%%100 ==0){
          print(rep)
        }
      }
      
      CI.BR['CP.LE'] = 0
      CI.BR['CP'] = 0
      CI.BR['CP.UE'] = 0
      CI.BR[which(CI.BR["lower"]<CI.BR["theta"] & CI.BR["upper"]>CI.BR["theta"]),"CP"]=1
      CI.BR[which(CI.BR["lower"]>CI.BR["theta"]),"CP.LE"]=1
      CI.BR[which(CI.BR["upper"]<CI.BR["theta"]),"CP.UE"]=1
      CI.BR["Length"] = CI.BR["upper"] - CI.BR["lower"]
      ACI.BR = colMeans(CI.BR[,c("CP.LE","CP","CP.UE","Length")])
      ACP.BR = as.numeric(ACI.BR["CP"])*100
      ACP.LE.BR = as.numeric(ACI.BR["CP.LE"])*100
      ACP.UE.BR = as.numeric(ACI.BR["CP.UE"])*100
      Ex.length.BR = as.numeric(ACI.BR["Length"])	
      
      str1 = paste("M=",M,"m=",m, "mu=",paste(muTotal,collapse = ' '),"k=",k, "nTotal=",paste(nTotal,collapse = ' '))
      str2 = paste("sigmaSqTotal=", paste(sigmaSqTotal,collapse = ' '),"deltaTotal=",paste(deltaTotal,collapse = ' '),"theta=",paste(theta,collapse = ' '))   
      
      str3 =paste( "LE|BR =",paste(ACP.LE.BR,collapse = ' '),"CP|BR =", paste(ACP.BR,collapse = ' '),"UE|BR =", paste(ACP.UE.BR,collapse = ' '), "AL|BR =", paste(Ex.length.BR,collapse = ' '))
      
      alloutput = paste(str1,str2,str3,collapse = ' ')
      folderName = unlist(strsplit(currentFolder, "//"))
      dir.create(file.path(resultoutput), showWarnings = FALSE)
      write(alloutput, file = paste(resultoutput,"/Results_BR.txt",sep = ""), append = TRUE)
    }


    ##############################################################################
    #                     Algorithm 3: GPQ                                       #
    ##############################################################################
    
    
    resultoutput = "./GPQ_result/"
    
    #Delete all result
    delfiles = dir(path=resultoutput ,pattern="*.txt")
    file.remove(file.path(resultoutput, delfiles))
    
    M = 5000
    m = 2500
    
    
    listFolder = list.dirs("./DataSource/", recursive=FALSE)
    
    eachFolder = listFolder[1]
    print(listFolder)
    for (eachFolder in listFolder){
      
      print(paste("Start Folder",eachFolder,sep=""))
      ## read data from config file
      currentFolder = eachFolder
      configfile = read.csv(file=paste(currentFolder,"/config.csv",sep = ""))
      nTotal = configfile$config_ntotal
      sigmaSqTotal = configfile$config_sigmaSq
      deltaTotal = configfile$config_delta
      muTotal = -sigmaSqTotal/2
      k = length(nTotal)
      
      
      # calculate parameter from data in config file
      deltanonTotal = 1 - deltaTotal
      size_nonzeroTotal = nTotal*(1-deltaTotal)
      Ex_LN = exp(muTotal+(sigmaSqTotal/2))
      cv_LN = sqrt(exp(sigmaSqTotal)-1)
      theta = log(deltanonTotal)+(2*muTotal+sigmaSqTotal)+log(exp(sigmaSqTotal)-deltanonTotal)
      theta12 = theta[1]-theta[2]
      theta13 = theta[1]-theta[3]
      theta23 = theta[2]-theta[3]
      
      #read summary data
      summaryData = read.csv(file= paste(currentFolder,"/basicSummary_K",k,".csv",sep = ""))
      
      #GPQ
      CI.GPQ = data.frame()
      D.GPQ = rep(0,m)
      
      rep = 1
      for(rep in seq(1:M)){
        
        #Select data from repeat loop number
        currentData = summaryData[which(summaryData$loop==rep),]
        
        #Extract each parameter 
        delta_hatnon = currentData$delta_hatnon
        delta_hat = 1-delta_hatnon
        size_nonzero = currentData$size_non
        size_zero = currentData$size_zero
        mu_hat = currentData$mu_hat
        sigmaSq_hat = currentData$sigmaSq_hat
        n = size_nonzero + size_zero
        
        thetahat = log(delta_hatnon) + (2*(mu_hat+sigmaSq_hat))
        thetahat.12 = thetahat[1]-thetahat[2]
        thetahat.13 = thetahat[1]-thetahat[3]
        thetahat.23 = thetahat[2]-thetahat[3]
        
        Var_varDLN1 = delta_hat/(n*delta_hatnon)
        Var_varDLN2 = (sigmaSq_hat/size_nonzero) + (2*sigmaSq_hat^2/(size_nonzero-1))
        Var_varDLN = Var_varDLN1 + (4*Var_varDLN2)
        
        Var_varDLN12 = Var_varDLN[1]+Var_varDLN[2]
        Var_varDLN13 = Var_varDLN[1]+Var_varDLN[3]
        Var_varDLN23 = Var_varDLN[2]+Var_varDLN[3]
        
        
        
       idx.GPQ = 1
        
        while(idx.GPQ <= m){
          
          U2_1 = rchisq(1,size_nonzero[1]-1)/(size_nonzero[1]-1)
          U2_2 = rchisq(1,size_nonzero[2]-1)/(size_nonzero[2]-1)
          U2_3 = rchisq(1,size_nonzero[3]-1)/(size_nonzero[3]-1)
          U2 = c(U2_1, U2_2, U2_3)
          
          mu_GPQ1 = mu_hat[1] - ((rnorm(1)/sqrt(U2[1]))*sqrt(sigmaSq_hat[1]/size_nonzero[1]))
          mu_GPQ2 = mu_hat[2] - ((rnorm(1)/sqrt(U2[2]))*sqrt(sigmaSq_hat[2]/size_nonzero[2]))
          mu_GPQ3 = mu_hat[3] - ((rnorm(1)/sqrt(U2[3]))*sqrt(sigmaSq_hat[3]/size_nonzero[3]))
          mu_GPQ = c(mu_GPQ1, mu_GPQ2, mu_GPQ3)
          sigmaSq_GPQ = sigmaSq_hat/U2
          
          delta_GPQ1 = (sin(asin(sqrt(delta_hat[1]))-(rnorm(1)/(2*sqrt(n[1])))))^2
          delta_GPQ2 = (sin(asin(sqrt(delta_hat[2]))-(rnorm(1)/(2*sqrt(n[2])))))^2
          delta_GPQ3 = (sin(asin(sqrt(delta_hat[3]))-(rnorm(1)/(2*sqrt(n[3])))))^2
          delta_GPQ4 = (sin(asin(sqrt(delta_hat[4]))-(rnorm(1)/(2*sqrt(n[4])))))^2
          delta_GPQ5 = (sin(asin(sqrt(delta_hat[5]))-(rnorm(1)/(2*sqrt(n[5])))))^2
          delta_GPQ = c(delta_GPQ1, delta_GPQ2, delta_GPQ3)
          
          theta_GPQ_i  = log(1-delta_GPQ) + 2*(mu_GPQ+ sigmaSq_GPQ)
          theta_GPQ12 = theta_GPQ_i[1]-theta_GPQ_i[2]
          theta_GPQ13 = theta_GPQ_i[1]-theta_GPQ_i[3]
          theta_GPQ23 = theta_GPQ_i[2]-theta_GPQ_i[3]
          
          D12.GPQ = abs((thetahat.12-theta_GPQ12)/sqrt(Var_varDLN12))
          D13.GPQ = abs((thetahat.13-theta_GPQ13)/sqrt(Var_varDLN13))
          D23.GPQ = abs((thetahat.23-theta_GPQ23)/sqrt(Var_varDLN23))
          
          D.GPQ[idx.GPQ] = max(D12.GPQ, D13.GPQ, D23.GPQ)
          idx.GPQ=idx.GPQ+1
        }
        
        q.GPQ =   quantile(D.GPQ,0.95,type=8)
        L.GPQ12 = thetahat.12-(q.GPQ*(sqrt(Var_varDLN12)))                         
        U.GPQ12 = thetahat.12+(q.GPQ*(sqrt(Var_varDLN12)))
        CI.GPQ12 = data.frame(rep,"CI.GPQ12",L.GPQ12,U.GPQ12,theta12)
        names(CI.GPQ12) = c("rep","CI","lower","upper","theta")
        CI.GPQ = rbind(CI.GPQ,CI.GPQ12)                                    
        
        L.GPQ13 = thetahat.13-(q.GPQ*(sqrt(Var_varDLN13)))                         
        U.GPQ13 = thetahat.13+(q.GPQ*(sqrt(Var_varDLN13)))
        CI.GPQ13 = data.frame(rep,"CI.GPQ13",L.GPQ13,U.GPQ13,theta13)
        names(CI.GPQ13) = c("rep","CI","lower","upper","theta")
        CI.GPQ = rbind(CI.GPQ,CI.GPQ13)                                        
        
        L.GPQ23 = thetahat.23-(q.GPQ*(sqrt(Var_varDLN23)))                         
        U.GPQ23 = thetahat.23+(q.GPQ*(sqrt(Var_varDLN23)))
        CI.GPQ23 = data.frame(rep,"CI.GPQ23",L.GPQ23,U.GPQ23,theta23)
        names(CI.GPQ23) = c("rep","CI","lower","upper","theta")
        CI.GPQ = rbind(CI.GPQ,CI.GPQ23)                                       
        
        
      }
      
      
      CI.GPQ['CP.LE'] = 0
      CI.GPQ['CP'] = 0
      CI.GPQ['CP.UE'] = 0
      CI.GPQ[which(CI.GPQ["lower"]<CI.GPQ["theta"] & CI.GPQ["upper"]>CI.GPQ["theta"]),"CP"]=1
      CI.GPQ[which(CI.GPQ["lower"]>CI.GPQ["theta"]),"CP.LE"]=1
      CI.GPQ[which(CI.GPQ["upper"]<CI.GPQ["theta"]),"CP.UE"]=1
      CI.GPQ["Length"] = CI.GPQ["upper"] - CI.GPQ["lower"]
      ACI.GPQ = colMeans(CI.GPQ[,c("CP.LE","CP","CP.UE","Length")])
      ACP.GPQ = as.numeric(ACI.GPQ["CP"])*100
      ACP.LE.GPQ = as.numeric(ACI.GPQ["CP.LE"])*100
      ACP.UE.GPQ = as.numeric(ACI.GPQ["CP.UE"])*100
      Ex.length.GPQ = as.numeric(ACI.GPQ["Length"])	  
      
      str1 = paste("M=",M,"m=",m, "mu=",paste(muTotal,collapse = ' '),"k=",k, "nTotal=",paste(nTotal,collapse = ' '))
      str2 = paste("sigmaSqTotal=", paste(sigmaSqTotal,collapse = ' '),"deltaTotal=",paste(deltaTotal,collapse = ' '),"theta=",paste(theta,collapse = ' '))   
      
      str3 =paste( "LE|GPQ =",paste(ACP.LE.GPQ,collapse = ' '),"CP|GPQ =", paste(ACP.GPQ,collapse = ' '),"UE|GPQ =", paste(ACP.UE.GPQ,collapse = ' '), "AL|GPQ =", paste(Ex.length.GPQ,collapse = ' '))
      
      alloutput = paste(str1,str2,str3,collapse = ' ')
      folderName = unlist(strsplit(currentFolder, "//"))
      dir.create(file.path(resultoutput), showWarnings = FALSE)
      write(alloutput, file = paste(resultoutput,"/Results_GPQ.txt",sep = ""), append = TRUE)
    }
    



    ##############################################################################
    #                  Algorithm 4: PB                                           #
    ##############################################################################    
    ### This function for loading the raw dataset
    loadXlist = function(folderselect,Kselect){
      D1 = read.csv(file=paste(folderselect,"/D",Kselect,"_xlist.csv",sep = ""),header = FALSE)
      return (D1)
    }
    
    
    resultoutput = "./PB_result/"
    
    #Delete all result
    delfiles = dir(path=resultoutput ,pattern="*.txt")
    file.remove(file.path(resultoutput, delfiles))
    
    M = 5000
    m = 2500
    
    
    listFolder = list.dirs("./DataSource/", recursive=FALSE)
    
    eachFolder = listFolder[1]
    print(listFolder)
    for (eachFolder in listFolder){
      
      print(paste("Start Folder",eachFolder,sep=""))
      ## read data from config file
      currentFolder = eachFolder
      configfile = read.csv(file=paste(currentFolder,"/config.csv",sep = ""))
      nTotal = configfile$config_ntotal
      sigmaSqTotal = configfile$config_sigmaSq
      deltaTotal = configfile$config_delta
      muTotal = -sigmaSqTotal/2
      k = length(nTotal)
      
      
      # calculate parameter from data in config file
      deltanonTotal = 1 - deltaTotal
      size_nonzeroTotal = nTotal*(1-deltaTotal)
      Ex_LN = exp(muTotal+(sigmaSqTotal/2))
      cv_LN = sqrt(exp(sigmaSqTotal)-1)
      theta = log(deltanonTotal)+(2*muTotal+sigmaSqTotal)+log(exp(sigmaSqTotal)-deltanonTotal)
      theta12 = theta[1]-theta[2]
      theta13 = theta[1]-theta[3]
      theta23 = theta[2]-theta[3]
      
      #read summary data
      summaryData = read.csv(file= paste(currentFolder,"/basicSummary_K",k,".csv",sep = ""))
      
      #PB: Parametric bootstrap
      CI.PB = data.frame()
      D.PB = rep(0,m)
      
      rep = 1
      dataframe.x1 = loadXlist(currentFolder,1)
      dataframe.x2 = loadXlist(currentFolder,2)
      dataframe.x3 = loadXlist(currentFolder,3)
      for(rep in seq(1:M)){
        
        #Select data from repeat loop number
        currentData = summaryData[which(summaryData$loop==rep),]
        
        #Extract each parameter 
        delta_hatnon = currentData$delta_hatnon
        delta_hat = 1-delta_hatnon
        size_nonzero = currentData$size_non
        size_zero = currentData$size_zero
        mu_hat = currentData$mu_hat
        sigmaSq_hat = currentData$sigmaSq_hat
        n = size_nonzero + size_zero
        
        thetahat = log(delta_hatnon) + (2*(mu_hat+sigmaSq_hat))
        thetahat.12 = thetahat[1]-thetahat[2]
        thetahat.13 = thetahat[1]-thetahat[3]
        thetahat.23 = thetahat[2]-thetahat[3]
        
        Var_varDLN1 = delta_hat/(n*delta_hatnon)
        Var_varDLN2 = (sigmaSq_hat/size_nonzero) + (2*sigmaSq_hat^2/(size_nonzero-1))
        Var_varDLN = Var_varDLN1 + (4*Var_varDLN2)
        
        Var_varDLN12 = Var_varDLN[1]+Var_varDLN[2]
        Var_varDLN13 = Var_varDLN[1]+Var_varDLN[3]
        Var_varDLN23 = Var_varDLN[2]+Var_varDLN[3]
        
        
        ##############################################################################    
        x1 = as.numeric(dataframe.x1[rep,])
        x2 = as.numeric(dataframe.x2[rep,])
        x3 = as.numeric(dataframe.x3[rep,])
        
        idx.PB = 1
        
        while(idx.PB <= m){
          
          #Bootstrap
          A1 = sample(1:n[1], size=n[1], replace=TRUE)
          X1.BS = x1[A1]
          nonzero_value1.BS = log(X1.BS[which(X1.BS!=0)])                          #find values is not zero
          size_nonzero1.BS =length(nonzero_value1.BS) 
          size_zero1.BS = n[1]-size_nonzero1.BS
          mu_hat1.BS = mean(nonzero_value1.BS)
          sigmaSq_hat1.BS = var(nonzero_value1.BS)
          delta_hat1.BS = 1-(size_nonzero1.BS/n[1])  
          
          A2 = sample(1:n[2], size=n[2], replace=TRUE)
          X2.BS = x2[A2]
          nonzero_value2.BS = log(X2.BS[which(X2.BS!=0)])                          #find values is not zero
          size_nonzero2.BS =length(nonzero_value2.BS) 
          size_zero2.BS = n[2]-size_nonzero2.BS
          mu_hat2.BS = mean(nonzero_value2.BS)
          sigmaSq_hat2.BS = var(nonzero_value2.BS)
          delta_hat2.BS = 1-(size_nonzero2.BS/n[2])  
          
          A3 = sample(1:n[3], size=n[3], replace=TRUE)
          X3.BS = x3[A3]
          nonzero_value3.BS = log(X3.BS[which(X3.BS!=0)])                          #find values is not zero
          size_nonzero3.BS =length(nonzero_value3.BS) 
          size_zero3.BS = n[3]-size_nonzero3.BS
          mu_hat3.BS = mean(nonzero_value3.BS)
          sigmaSq_hat3.BS = var(nonzero_value3.BS)
          delta_hat3.BS = 1-(size_nonzero3.BS/n[3])
          
          
          size_nonzero.BS = c(size_nonzero1.BS, size_nonzero2.BS, size_nonzero3.BS)
          mu_hat.BS = c(mu_hat1.BS, mu_hat2.BS, mu_hat3.BS)
          sigmaSq_hat.BS = c(sigmaSq_hat1.BS, sigmaSq_hat2.BS, sigmaSq_hat3.BS)
          delta_hat.BS = c(delta_hat1.BS, delta_hat2.BS, delta_hat3.BS)
          delta_hatnon.BS = 1 - delta_hat.BS
          thetahat.BS = log(delta_hatnon.BS) + (2*(mu_hat.BS+sigmaSq_hat.BS))
          
          
          #Parametric boostrap
          mu_hat1.PB =  mu_hat1.BS + (rnorm(1)*sqrt(sigmaSq_hat1.BS/size_nonzero1.BS))
          mu_hat2.PB =  mu_hat2.BS + (rnorm(1)*sqrt(sigmaSq_hat2.BS/size_nonzero2.BS))
          mu_hat3.PB =  mu_hat3.BS + (rnorm(1)*sqrt(sigmaSq_hat3.BS/size_nonzero3.BS))
          
          sigmaSq_hat1.PB = sigmaSq_hat1.BS*rchisq(1,size_nonzero1.BS-1, ncp = 0)/(size_nonzero1.BS-1)
          sigmaSq_hat2.PB = sigmaSq_hat2.BS*rchisq(1,size_nonzero2.BS-1, ncp = 0)/(size_nonzero2.BS-1)
          sigmaSq_hat3.PB = sigmaSq_hat3.BS*rchisq(1,size_nonzero3.BS-1, ncp = 0)/(size_nonzero3.BS-1)
          
          delta_hat1.PB = rbeta(1, size_zero1.BS+0.5, size_nonzero1.BS+0.5)
          delta_hat2.PB = rbeta(1, size_zero2.BS+0.5, size_nonzero2.BS+0.5)
          delta_hat3.PB = rbeta(1, size_zero3.BS+0.5, size_nonzero3.BS+0.5)
          delta_hat.PB = c(delta_hat1.PB, delta_hat2.PB, delta_hat3.PB)
          delta_hatnon.PB = 1-delta_hat.PB
          size_nonzero.PB = n*delta_hatnon.PB
          mu_hat.PB = c(mu_hat1.PB, mu_hat2.PB, mu_hat3.PB)
          sigmaSq_hat.PB = c(sigmaSq_hat1.PB, sigmaSq_hat2.PB, sigmaSq_hat3.PB)
          thetahat.PB = log(delta_hatnon.PB) + (2*(mu_hat.PB+sigmaSq_hat.PB))
          
          Var_varDLN1.PB = delta_hat.PB/(n*delta_hatnon.PB)
          Var_varDLN2.PB = (sigmaSq_hat.PB/size_nonzero.PB) + (2*sigmaSq_hat.PB^2/(size_nonzero.PB-1))
          Var_varDLN.PB = Var_varDLN1.PB + (4*Var_varDLN2.PB) 
          V12.PB = Var_varDLN.PB[1] + Var_varDLN.PB[2]
          V13.PB = Var_varDLN.PB[1] + Var_varDLN.PB[3]
          V23.PB = Var_varDLN.PB[2] + Var_varDLN.PB[3]
          
          
          D12.PB  = abs(((thetahat.PB[1] - thetahat.PB[2])-(thetahat.BS[1] - thetahat.BS[2]))/sqrt(V12.PB))
          D13.PB  = abs(((thetahat.PB[1] - thetahat.PB[3])-(thetahat.BS[1] - thetahat.BS[3]))/sqrt(V13.PB))
          D23.PB  = abs(((thetahat.PB[2] - thetahat.PB[3])-(thetahat.BS[2] - thetahat.BS[3]))/sqrt(V23.PB))
          
          D.PB[idx.PB] = max(D12.PB,D13.PB,D23.PB)
          if(is.na(D.PB[idx.PB])){
            next
          }
          
          idx.PB=idx.PB+1
          
          
        }
        
        q.PB = quantile(D.PB,0.95,type=8)
        L.PB12 = thetahat.12-(q.PB*(sqrt(Var_varDLN12)))                        
        U.PB12 = thetahat.12+(q.PB*(sqrt(Var_varDLN12)))  
        CI.PB12 = data.frame(rep,"CI.PB12",L.PB12,U.PB12,theta12)
        names(CI.PB12) = c("rep","CI","lower","upper","theta")
        CI.PB = rbind(CI.PB,CI.PB12)
        
        L.PB13 = thetahat.13-(q.PB*(sqrt(Var_varDLN13)))                         
        U.PB13 = thetahat.13+(q.PB*(sqrt(Var_varDLN13)))                        
        CI.PB13 = data.frame(rep,"CI.PB13",L.PB13,U.PB13,theta13)
        names(CI.PB13) = c("rep","CI","lower","upper","theta")
        CI.PB = rbind(CI.PB,CI.PB13)                                   
        
        L.PB23 = thetahat.23-(q.PB*(sqrt(Var_varDLN23)))                         
        U.PB23 = thetahat.23+(q.PB*(sqrt(Var_varDLN23)))                         
        CI.PB23 = data.frame(rep,"CI.PB23",L.PB23,U.PB23,theta23)
        names(CI.PB23) = c("rep","CI","lower","upper","theta")
        CI.PB = rbind(CI.PB,CI.PB23)                                    
        
        
        
      }
      
      CI.PB['CP.LE'] = 0
      CI.PB['CP'] = 0
      CI.PB['CP.UE'] = 0
      CI.PB[which(CI.PB["lower"]<CI.PB["theta"] & CI.PB["upper"]>CI.PB["theta"]),"CP"]=1
      CI.PB[which(CI.PB["lower"]>CI.PB["theta"]),"CP.LE"]=1
      CI.PB[which(CI.PB["upper"]<CI.PB["theta"]),"CP.UE"]=1
      CI.PB["Length"] = CI.PB["upper"] - CI.PB["lower"]
      ACI.PB = colMeans(CI.PB[,c("CP.LE","CP","CP.UE","Length")])
      ACP.PB = as.numeric(ACI.PB["CP"])*100
      ACP.LE.PB = as.numeric(ACI.PB["CP.LE"])*100
      ACP.UE.PB = as.numeric(ACI.PB["CP.UE"])*100
      Ex.length.PB = as.numeric(ACI.PB["Length"])	
      
      str1 = paste("M=",M,"m=",m, "mu=",paste(muTotal,collapse = ' '),"k=",k, "nTotal=",paste(nTotal,collapse = ' '))
      str2 = paste("sigmaSqTotal=", paste(sigmaSqTotal,collapse = ' '),"deltaTotal=",paste(deltaTotal,collapse = ' '),"theta=",paste(theta,collapse = ' '))   
      
      str3 =paste( "LE|PB =",paste(ACP.LE.PB,collapse = ' '),"CP|PB =", paste(ACP.PB,collapse = ' '),"UE|PB =", paste(ACP.UE.PB,collapse = ' '), "AL|PB =", paste(Ex.length.PB,collapse = ' '))
      
      alloutput = paste(str1,str2,str3,collapse = ' ')
      folderName = unlist(strsplit(currentFolder, "//"))
      dir.create(file.path(resultoutput), showWarnings = FALSE)
      write(alloutput, file = paste(resultoutput,"/Results_PB.txt",sep = ""), append = TRUE)
    }
    



    ##################################################################################
    ##################################################################################
    # Five sample groups
    # k = 5
    ##################################################################################
    ##################################################################################
    
      #Simulation runs
      M = 5000
      m = 2500
      
      # Function for generating random sample drawn from delta-lognormal distributions with k-th populations
      library(EnvStats)
      
      xGen = function(k,nTotal,Ex_LN,cv_LN,deltaTotal){
        sum_nonzero_value2 = rep(0,k)
        mu_hat = rep(0,k)
        sigmaSq_hat= rep(0,k)
        delta_hatnon= rep(0,k)
        size_non = rep(0,k)
        size_zero = rep(0,k)
        xlist = list()
        
        
        for (l in 1:k){
          #print(l)
          while(TRUE){
            x = rzmlnormAlt(nTotal[l], Ex_LN[l], cv_LN[l], deltaTotal[l])
            xlist[l] = list(x)
            
            Temp.nonzero_value = log(x[which(x!=0)])                         
            Temp.size_nonzero =length(Temp.nonzero_value)                          
            if(Temp.size_nonzero < (nTotal[l]-nTotal[l]*deltaTotal[l])){
              next
            }
            nonzero_value = Temp.nonzero_value
            size_nonzero = Temp.size_nonzero
            sum_nonzero_value2[l] = sum(nonzero_value^2)
            mu_hat[l] = mean(nonzero_value)
            sigmaSq_hat[l] = var(nonzero_value)
            
            delta_hatnon[l] = size_nonzero/nTotal[l]
            size_non[l] = size_nonzero
            size_zero[l] = nTotal[l]-size_nonzero
            
            break
          }
        }
        return(list("data" = cbind(size_non, size_zero, mu_hat, sigmaSq_hat, delta_hatnon, sum_nonzero_value2),"xlist"=xlist))
      }
    
    
    
      # Parameter estimations
      writeDataToFile=function(x,k,rep,currentFolderPath){
        datatable = t(data.frame(x))
        filename = paste(currentFolderPath,'/D',k,'_xlist.csv',sep="")
        write.table( datatable, filename  , append= T, sep=',',col.names=FALSE,row.names = FALSE )
      }
      
      
      createDatasource = function(nCase,sCase,dCase){
      
      ###############################Delta Varied########################################
      folderSourceName = paste("n",nCase,"s",sCase,"d",dCase, sep = "")
      print(folderSourceName)
      ###################################################################################
      #Create folder
      dir.create(file.path("./DataSource/"), showWarnings = FALSE)
      
      ## set current folder to keep all data  ./Datasource/Data_m_2500_n1_100n2_100.......
      currentFolderPath = paste("./DataSource/",folderSourceName,sep="")
      dir.create(file.path(currentFolderPath),showWarnings = FALSE)
      
      return (currentFolderPath)
      
      }
    
    
      ntotal_allcase = matrix(c( 30,   50,   50,  100,  200,
                                 50,   50,   50,  100,  100,
                                 70,  100,  100,  200,  200), ncol = 5,byrow = TRUE)
      
      sigmaSq_allcase = matrix(c( 1,  1,    2,  2,  3), ncol = 5,byrow = TRUE) 
      
      delta_allcase = matrix(c( 0.1,  0.1,  0.2,  0.2,  0.2,
                                0.2,  0.2,  0.3,  0.3,  0.5,
                                0.2,  0.3,  0.5,  0.5,  0.7,
                                0.5,  0.5,  0.5,  0.7,  0.7), ncol = 5,byrow = TRUE) 
    
    
    for (id_ntotal_alcase in seq(1:nrow(ntotal_allcase))){
      config_ntotal =ntotal_allcase[id_ntotal_alcase,]
      
      for (id_sigmaSq_allcase in seq(1:nrow(sigmaSq_allcase))){
        config_sigmaSq =sigmaSq_allcase[id_sigmaSq_allcase,]
        
        for (id_delta_allcase in seq(1:nrow(delta_allcase))){
          config_delta =delta_allcase[id_delta_allcase,]
          
          currentFolderPath = createDatasource(id_ntotal_alcase,id_sigmaSq_allcase,id_delta_allcase)
          
          allData = data.frame(data.frame(config_ntotal),data.frame(config_sigmaSq),data.frame(config_delta))
          #data_config_object = data.frame(M, m, n1, n2, mu1, mu2, sigmaSq1, sigmaSq2, delta1, delta2, alpha)
          write.csv(allData,file=paste(currentFolderPath,"/config.csv",sep=""),row.names = FALSE)
        }
        
      }
      
    }
    
    listFolder = list.dirs("./DataSource/", recursive=FALSE)
    
    eachFolder = listFolder[1]
    print(listFolder)
    for (eachFolder in listFolder){
      currentFolder = eachFolder
      configfile = read.csv(file=paste(currentFolder,"/config.csv",sep = ""))
      nTotal = configfile$config_ntotal
      sigmaSqTotal = configfile$config_sigmaSq
      deltaTotal = configfile$config_delta
      
      
      k = length(nTotal)
      
      muTotal =  -sigmaSqTotal/2
      mu1 = muTotal[1]
      mu2 = muTotal[2]
      mu3 = muTotal[3]
      mu4 = muTotal[4]
      mu5 = muTotal[5]
      
      deltanonTotal = 1 - deltaTotal
      size_nonzeroTotal = nTotal*(1-deltaTotal)
      
      
      Ex_LN = exp(muTotal+(sigmaSqTotal/2))
      cv_LN = sqrt(exp(sigmaSqTotal)-1)
      theta = log(deltanonTotal)+(2*muTotal+sigmaSqTotal)+log(exp(sigmaSqTotal)-deltanonTotal)
      theta12 = theta[1]-theta[2]
      theta13 = theta[1]-theta[3]
      theta14 = theta[1]-theta[4]
      theta15 = theta[1]-theta[5]
      theta23 = theta[2]-theta[3]
      theta24 = theta[2]-theta[4]
      theta25 = theta[2]-theta[5]
      theta34 = theta[3]-theta[4]
      theta35 = theta[3]-theta[5]
      theta45 = theta[4]-theta[5]
      
      
      
      AllDataset = data.frame()
      ## check file .xlxs is exist it will be deleted
      
      delfiles = dir(path=currentFolder ,pattern="D.*xlist.csv")
      #Check its existence
      file.remove(file.path(currentFolder, delfiles))
      
      idx_rep =1
      for ( idx_rep in seq(1:M)){
        #Start generate data
        packData = xGen(k, nTotal, Ex_LN, cv_LN, deltaTotal)
        dataset= data.frame(loop = rep(idx_rep,k),num_k= 1:k)
        data  = packData$data
        xlist = packData$xlist
        size_nonzero = data
        size_zero = data[,2]
        mu_hat = data[,3]
        sigmaSq_hat  = data[,4]
        delta_hatnon = data[,5]
        delta_hat  = 1-delta_hatnon
        sum_nonzero_value2 = data[,6]
        dataset = cbind(dataset, data)
        
        ## store Dataset each loop
        AllDataset = rbind(AllDataset,dataset)
        
        
        # write raw xlist data
        for (idx_xlist in seq(1:length(xlist))){
          eachlist = xlist[idx_xlist]
          writeDataToFile(eachlist,idx_xlist,idx_rep,currentFolder)
        }
        
      }
      write.table( AllDataset, paste(currentFolder,"/basicSummary_K",k,".csv",sep="")  , sep=',',row.names = FALSE )
    }
    
    
    ##############################################################################
    #                 Algorithm 1: B-PMB                                         #
    ##############################################################################  
    library(rstan)            #Used to read stan file
    
    ### This function for loading the raw dataset
    loadXlist = function(folderselect,Kselect){
      D1 = read.csv(file=paste(folderselect,"/D",Kselect,"_xlist.csv",sep = ""),header = FALSE)
      return (D1)
    }
    
    #rstan_options(auto_write=TRUE)
    
    resultoutput = "./BPM_result/"
    
    #Delete all result
    delfiles = dir(path=resultoutput ,pattern="*.txt")
    file.remove(file.path(resultoutput, delfiles))
    
    M = 5000
    m = 2500
    
    
    listFolder = list.dirs("./DataSource/", recursive=FALSE)
    
    eachFolder = listFolder[1]
    print(listFolder)
    for (eachFolder in listFolder){
      
      print(paste("Start Folder",eachFolder,sep=""))
      ## read data from config file
      currentFolder = eachFolder
      configfile = read.csv(file=paste(currentFolder,"/config.csv",sep = ""))
      nTotal = configfile$config_ntotal
      sigmaSqTotal = configfile$config_sigmaSq
      deltaTotal = configfile$config_delta
      muTotal = -sigmaSqTotal/2
      k = length(nTotal)
      
      
      # calculate parameter from data in config file
      deltanonTotal = 1 - deltaTotal
      size_nonzeroTotal = nTotal*(1-deltaTotal)
      Ex_LN = exp(muTotal+(sigmaSqTotal/2))
      cv_LN = sqrt(exp(sigmaSqTotal)-1)
      theta = log(deltanonTotal)+(2*muTotal+sigmaSqTotal)+log(exp(sigmaSqTotal)-deltanonTotal)
      theta12 = theta[1]-theta[2]
      theta13 = theta[1]-theta[3]
      theta14 = theta[1]-theta[4]
      theta15 = theta[1]-theta[5]
      theta23 = theta[2]-theta[3]
      theta24 = theta[2]-theta[4]
      theta25 = theta[2]-theta[5]
      theta34 = theta[3]-theta[4]
      theta35 = theta[3]-theta[5]
      theta45 = theta[4]-theta[5]
      
      #read summary data
      summaryData = read.csv(file= paste(currentFolder,"/basicSummary_K",k,".csv",sep = ""))
      
      #BPM: Bayesian-based probability-matching prior 
      CI.BPM = data.frame()
      D.BPM  = rep(0,m)
      
      rep = 1
      dataframe.x1 = loadXlist(currentFolder,1)
      dataframe.x2 = loadXlist(currentFolder,2)
      dataframe.x3 = loadXlist(currentFolder,3)
      dataframe.x4 = loadXlist(currentFolder,4)
      dataframe.x5 = loadXlist(currentFolder,5)
      
      for(rep in seq(1:M)){
        
        #Select data from repeat loop number
        currentData = summaryData[which(summaryData$loop==rep),]
        
        #Extract each parameter 
        delta_hatnon = currentData$delta_hatnon
        delta_hat = 1-delta_hatnon
        size_nonzero = currentData$size_non
        size_zero = currentData$size_zero
        mu_hat = currentData$mu_hat
        sigmaSq_hat = currentData$sigmaSq_hat
        n = size_nonzero + size_zero
        
        thetahat = log(delta_hatnon) + (2*(mu_hat+sigmaSq_hat))
        thetahat.12 = thetahat[1]-thetahat[2]
        thetahat.13 = thetahat[1]-thetahat[3]
        thetahat.14 = thetahat[1]-thetahat[4]
        thetahat.15 = thetahat[1]-thetahat[5]
        thetahat.23 = thetahat[2]-thetahat[3]
        thetahat.24 = thetahat[2]-thetahat[4]
        thetahat.25 = thetahat[2]-thetahat[5]
        thetahat.34 = thetahat[3]-thetahat[4]
        thetahat.35 = thetahat[3]-thetahat[5]
        thetahat.45 = thetahat[4]-thetahat[5]
        
        Var_varDLN1 = delta_hat/(n*delta_hatnon)
        Var_varDLN2 = (sigmaSq_hat/size_nonzero) + (2*sigmaSq_hat^2/(size_nonzero-1))
        Var_varDLN = Var_varDLN1 + (4*Var_varDLN2)
        
        Var_varDLN12 = Var_varDLN[1]+Var_varDLN[2]
        Var_varDLN13 = Var_varDLN[1]+Var_varDLN[3]
        Var_varDLN14 = Var_varDLN[1]+Var_varDLN[4]
        Var_varDLN15 = Var_varDLN[1]+Var_varDLN[5]
        Var_varDLN23 = Var_varDLN[2]+Var_varDLN[3]
        Var_varDLN24 = Var_varDLN[2]+Var_varDLN[4]
        Var_varDLN25 = Var_varDLN[2]+Var_varDLN[5]
        Var_varDLN34 = Var_varDLN[3]+Var_varDLN[4]
        Var_varDLN35 = Var_varDLN[3]+Var_varDLN[5]
        Var_varDLN45 = Var_varDLN[4]+Var_varDLN[5]
        
         
        
        ############################### BPM ######################################## 
        ## chain = 4, iter = 1250 to obtain m = 4*(1250/2) = 2500  
        ## refresh = 0 denoted as Stan running without output details
        
        x1 = as.numeric(dataframe.x1[rep,])
        x2 = as.numeric(dataframe.x2[rep,])
        x3 = as.numeric(dataframe.x3[rep,])
        x4 = as.numeric(dataframe.x4[rep,])
        x5 = as.numeric(dataframe.x5[rep,])
        
        
        fit.pm = stan('Stan_PMBprior_K5.stan', data=list(n1=n[1], sizenonzero1=size_nonzero[1], 
                                                                 y1=x1, sizezero1=size_zero[1], mu1=mu_hat[1], 
                                                                 sigmaSq1=sigmaSq_hat[1], delta1=delta_hat[1],
                                                                 
                                                                 n2=n[2], sizenonzero2=size_nonzero[2], 
                                                                 y2=x2, sizezero2=size_zero[2], mu2=mu_hat[2], 
                                                                 sigmaSq2=sigmaSq_hat[2], delta2=delta_hat[2],
                                                                 
                                                                 n3=n[3], sizenonzero3=size_nonzero[3], 
                                                                 y3=x3, sizezero3=size_zero[3], mu3=mu_hat[3], 
                                                                 sigmaSq3=sigmaSq_hat[3], delta3=delta_hat[3],
                                                                 
                                                                 n4=n[4], sizenonzero4=size_nonzero[4], 
                                                                 y4=x4, sizezero4=size_zero[4], mu4=mu_hat[4], 
                                                                 sigmaSq4=sigmaSq_hat[4], delta4=delta_hat[4], 
                                                                 
                                                                 n5=n[5], sizenonzero5=size_nonzero[5], 
                                                                 y5=x5, sizezero5=size_zero[5], mu5=mu_hat[5], 
                                                                 sigmaSq5=sigmaSq_hat[5], delta5=delta_hat[5]), 
                      iter = 1250, refresh = 0)
        
        # Extracting the posterior draws
        post.mu.pm1 = extract(fit.pm)$mu1
        post.sigmaSq.pm1 = extract(fit.pm)$sigmaSq1
        post.delta.pm1 = extract(fit.pm)$delta1
        post_varDLN.pm1 = log(1-post.delta.pm1)+2*(post.mu.pm1+post.sigmaSq.pm1)
        
        post.mu.pm2 = extract(fit.pm)$mu2
        post.sigmaSq.pm2 = extract(fit.pm)$sigmaSq2
        post.delta.pm2 = extract(fit.pm)$delta2
        post_varDLN.pm2 = log(1-post.delta.pm2)+2*(post.mu.pm2+post.sigmaSq.pm2)
        
        post.mu.pm3 = extract(fit.pm)$mu3
        post.sigmaSq.pm3 = extract(fit.pm)$sigmaSq3
        post.delta.pm3 = extract(fit.pm)$delta3
        post_varDLN.pm3 = log(1-post.delta.pm3)+2*(post.mu.pm3+post.sigmaSq.pm3)
        
        post.mu.pm4 = extract(fit.pm)$mu4
        post.sigmaSq.pm4 = extract(fit.pm)$sigmaSq4
        post.delta.pm4 = extract(fit.pm)$delta4
        post_varDLN.pm4 = log(1-post.delta.pm4)+2*(post.mu.pm4+post.sigmaSq.pm4)
        
        post.mu.pm5 = extract(fit.pm)$mu5
        post.sigmaSq.pm5 = extract(fit.pm)$sigmaSq5
        post.delta.pm5 = extract(fit.pm)$delta5
        post_varDLN.pm5 = log(1-post.delta.pm5)+2*(post.mu.pm5+post.sigmaSq.pm5)
        
        post_varDLN.pm12 = post_varDLN.pm1-post_varDLN.pm2
        post_varDLN.pm13 = post_varDLN.pm1-post_varDLN.pm3
        post_varDLN.pm14 = post_varDLN.pm1-post_varDLN.pm4
        post_varDLN.pm15 = post_varDLN.pm1-post_varDLN.pm5
        post_varDLN.pm23 = post_varDLN.pm2-post_varDLN.pm3
        post_varDLN.pm24 = post_varDLN.pm2-post_varDLN.pm4
        post_varDLN.pm25 = post_varDLN.pm2-post_varDLN.pm5
        post_varDLN.pm34 = post_varDLN.pm3-post_varDLN.pm4
        post_varDLN.pm35 = post_varDLN.pm3-post_varDLN.pm5
        post_varDLN.pm45 = post_varDLN.pm4-post_varDLN.pm5
        
        for(idx.BPM in 1:m){
          D12.BPM = abs(thetahat.12 - post_varDLN.pm12[idx.BPM])
          D13.BPM = abs(thetahat.13 - post_varDLN.pm13[idx.BPM])
          D14.BPM = abs(thetahat.14 - post_varDLN.pm14[idx.BPM])
          D15.BPM = abs(thetahat.15 - post_varDLN.pm15[idx.BPM])
          D23.BPM = abs(thetahat.23 - post_varDLN.pm23[idx.BPM])
          D24.BPM = abs(thetahat.24 - post_varDLN.pm24[idx.BPM])
          D25.BPM = abs(thetahat.25 - post_varDLN.pm25[idx.BPM])
          D34.BPM = abs(thetahat.34 - post_varDLN.pm34[idx.BPM])
          D35.BPM = abs(thetahat.35 - post_varDLN.pm35[idx.BPM])
          D45.BPM = abs(thetahat.45 - post_varDLN.pm45[idx.BPM])
          
          A.BPM = max(D12.BPM, D13.BPM, D14.BPM, D15.BPM, D23.BPM, D24.BPM, D25.BPM, D34.BPM, D35.BPM, D45.BPM)
          B.BPM = min(D12.BPM, D13.BPM, D14.BPM, D15.BPM, D23.BPM, D24.BPM, D25.BPM, D34.BPM, D35.BPM, D45.BPM)
          
          D.BPM[idx.BPM] = A.BPM - B.BPM
          
        }
        
        
        q.BPM =   quantile(D.BPM,0.95,type=8)
        L.BPM12 = thetahat.12 - q.BPM                     
        U.BPM12 = thetahat.12 + q.BPM
        CI.BPM12 = data.frame(rep,"CI.BPM12",L.BPM12,U.BPM12,theta12)
        names(CI.BPM12) = c("rep","CI","lower","upper","theta")
        CI.BPM = rbind(CI.BPM,CI.BPM12) 
        
        L.BPM13 = thetahat.13 - q.BPM                       
        U.BPM13 = thetahat.13 + q.BPM
        CI.BPM13 = data.frame(rep,"CI.BPM13",L.BPM13,U.BPM13,theta13)
        names(CI.BPM13) = c("rep","CI","lower","upper","theta")
        CI.BPM = rbind(CI.BPM,CI.BPM13) 
        
        L.BPM14 = thetahat.14 - q.BPM                       
        U.BPM14 = thetahat.14 + q.BPM
        CI.BPM14 = data.frame(rep,"CI.BPM14",L.BPM14,U.BPM14,theta14)
        names(CI.BPM14) = c("rep","CI","lower","upper","theta")
        CI.BPM = rbind(CI.BPM,CI.BPM14) 
        
        L.BPM15 = thetahat.15 - q.BPM                        
        U.BPM15 = thetahat.15 + q.BPM
        CI.BPM15 = data.frame(rep,"CI.BPM15",L.BPM15,U.BPM15,theta15)
        names(CI.BPM15) = c("rep","CI","lower","upper","theta")
        CI.BPM = rbind(CI.BPM,CI.BPM15) 
        
        L.BPM23 = thetahat.23 - q.BPM                    
        U.BPM23 = thetahat.23 + q.BPM
        CI.BPM23 = data.frame(rep,"CI.BPM23",L.BPM23,U.BPM23,theta23)
        names(CI.BPM23) = c("rep","CI","lower","upper","theta")
        CI.BPM = rbind(CI.BPM,CI.BPM23) 
        
        L.BPM24 = thetahat.24 - q.BPM                       
        U.BPM24 = thetahat.24 + q.BPM
        CI.BPM24 = data.frame(rep,"CI.BPM24",L.BPM24,U.BPM24,theta24)
        names(CI.BPM24) = c("rep","CI","lower","upper","theta")
        CI.BPM = rbind(CI.BPM,CI.BPM24) 
        
        L.BPM25 = thetahat.25 - q.BPM                        
        U.BPM25 = thetahat.25 + q.BPM
        CI.BPM25 = data.frame(rep,"CI.BPM25",L.BPM25,U.BPM25,theta25)
        names(CI.BPM25) = c("rep","CI","lower","upper","theta")
        CI.BPM = rbind(CI.BPM,CI.BPM25) 
        
        L.BPM34 = thetahat.34 - q.BPM                        
        U.BPM34 = thetahat.34 + q.BPM
        CI.BPM34 = data.frame(rep,"CI.BPM34",L.BPM34,U.BPM34,theta34)
        names(CI.BPM34) = c("rep","CI","lower","upper","theta")
        CI.BPM = rbind(CI.BPM,CI.BPM34) 
        
        L.BPM35 = thetahat.35 - q.BPM                        
        U.BPM35 = thetahat.35 + q.BPM
        CI.BPM35 = data.frame(rep,"CI.BPM35",L.BPM35,U.BPM35,theta35)
        names(CI.BPM35) = c("rep","CI","lower","upper","theta")
        CI.BPM = rbind(CI.BPM,CI.BPM35) 
        
        L.BPM45 = thetahat.45 - q.BPM                        
        U.BPM45 = thetahat.45 + q.BPM
        CI.BPM45 = data.frame(rep,"CI.BPM45",L.BPM45,U.BPM45,theta45)
        names(CI.BPM45) = c("rep","CI","lower","upper","theta")
        CI.BPM = rbind(CI.BPM,CI.BPM45)
        
        if(rep%%100 ==0){
          print(rep)
        }
      }
      
      CI.BPM['CP.LE'] = 0
      CI.BPM['CP'] = 0
      CI.BPM['CP.UE'] = 0
      CI.BPM[which(CI.BPM["lower"]<CI.BPM["theta"] & CI.BPM["upper"]>CI.BPM["theta"]),"CP"]=1
      CI.BPM[which(CI.BPM["lower"]>CI.BPM["theta"]),"CP.LE"]=1
      CI.BPM[which(CI.BPM["upper"]<CI.BPM["theta"]),"CP.UE"]=1
      CI.BPM["Length"] = CI.BPM["upper"] - CI.BPM["lower"]
      ACI.BPM = colMeans(CI.BPM[,c("CP.LE","CP","CP.UE","Length")])
      ACP.BPM = as.numeric(ACI.BPM["CP"])*100
      ACP.LE.BPM = as.numeric(ACI.BPM["CP.LE"])*100
      ACP.UE.BPM = as.numeric(ACI.BPM["CP.UE"])*100
      Ex.length.BPM = as.numeric(ACI.BPM["Length"])	
      
      str1 = paste("M=",M,"m=",m, "mu=",paste(muTotal,collapse = ' '),"k=",k, "nTotal=",paste(nTotal,collapse = ' '))
      str2 = paste("sigmaSqTotal=", paste(sigmaSqTotal,collapse = ' '),"deltaTotal=",paste(deltaTotal,collapse = ' '),"theta=",paste(theta,collapse = ' '))   
      
      str3 =paste( "LE|BPM =",paste(ACP.LE.BPM,collapse = ' '),"CP|BPM =", paste(ACP.BPM,collapse = ' '),"UE|BPM =", paste(ACP.UE.BPM,collapse = ' '), "AL|BPM =", paste(Ex.length.BPM,collapse = ' '))
      
      alloutput = paste(str1,str2,str3,collapse = ' ')
      folderName = unlist(strsplit(currentFolder, "//"))
      dir.create(file.path(resultoutput), showWarnings = FALSE)
      write(alloutput, file = paste(resultoutput,"/Results_BPM.txt",sep = ""), append = TRUE)
    }
    
    

    
    ##############################################################################
    #                 Algorithm 2 : B-RB                                         #
    ############################################################################## 
    library(rstan)            #Used to read stan file
    
    ### This function for loading the raw dataset
    loadXlist = function(folderselect,Kselect){
      D1 = read.csv(file=paste(folderselect,"/D",Kselect,"_xlist.csv",sep = ""),header = FALSE)
      return (D1)
    }
    
    resultoutput = "./BR_result/"
    
    #Delete all result
    delfiles = dir(path=resultoutput ,pattern="*.txt")
    file.remove(file.path(resultoutput, delfiles))
    
    M = 5000
    m = 2500
    
    
    listFolder = list.dirs("./DataSource/", recursive=FALSE)
    
    eachFolder = listFolder[1]
    print(listFolder)
    for (eachFolder in listFolder){
      
      print(paste("Start Folder",eachFolder,sep=""))
      ## read data from config file
      currentFolder = eachFolder
      configfile = read.csv(file=paste(currentFolder,"/config.csv",sep = ""))
      nTotal = configfile$config_ntotal
      sigmaSqTotal = configfile$config_sigmaSq
      deltaTotal = configfile$config_delta
      muTotal = -sigmaSqTotal/2
      k = length(nTotal)
      
      
      # calculate parameter from data in config file
      deltanonTotal = 1 - deltaTotal
      size_nonzeroTotal = nTotal*(1-deltaTotal)
      Ex_LN = exp(muTotal+(sigmaSqTotal/2))
      cv_LN = sqrt(exp(sigmaSqTotal)-1)
      theta = log(deltanonTotal)+(2*muTotal+sigmaSqTotal)+log(exp(sigmaSqTotal)-deltanonTotal)
      theta12 = theta[1]-theta[2]
      theta13 = theta[1]-theta[3]
      theta14 = theta[1]-theta[4]
      theta15 = theta[1]-theta[5]
      theta23 = theta[2]-theta[3]
      theta24 = theta[2]-theta[4]
      theta25 = theta[2]-theta[5]
      theta34 = theta[3]-theta[4]
      theta35 = theta[3]-theta[5]
      theta45 = theta[4]-theta[5]
      
      #read summary data
      summaryData = read.csv(file= paste(currentFolder,"/basicSummary_K",k,".csv",sep = ""))
      
      #BR: Bayesian-based reference prior 
      CI.BR = data.frame()
      D.BR  = rep(0,m)
      
      rep = 1
      dataframe.x1 = loadXlist(currentFolder,1)
      dataframe.x2 = loadXlist(currentFolder,2)
      dataframe.x3 = loadXlist(currentFolder,3)
      dataframe.x4 = loadXlist(currentFolder,4)
      dataframe.x5 = loadXlist(currentFolder,5)
      for(rep in seq(1:M)){
        
        #Select data from repeat loop number
        currentData = summaryData[which(summaryData$loop==rep),]
        
        #Extract each parameter 
        delta_hatnon = currentData$delta_hatnon
        delta_hat = 1-delta_hatnon
        size_nonzero = currentData$size_non
        size_zero = currentData$size_zero
        mu_hat = currentData$mu_hat
        sigmaSq_hat = currentData$sigmaSq_hat
        n = size_nonzero + size_zero
        
        thetahat = log(delta_hatnon) + (2*(mu_hat+sigmaSq_hat))
        thetahat.12 = thetahat[1]-thetahat[2]
        thetahat.13 = thetahat[1]-thetahat[3]
        thetahat.14 = thetahat[1]-thetahat[4]
        thetahat.15 = thetahat[1]-thetahat[5]
        thetahat.23 = thetahat[2]-thetahat[3]
        thetahat.24 = thetahat[2]-thetahat[4]
        thetahat.25 = thetahat[2]-thetahat[5]
        thetahat.34 = thetahat[3]-thetahat[4]
        thetahat.35 = thetahat[3]-thetahat[5]
        thetahat.45 = thetahat[4]-thetahat[5]
        
        Var_varDLN1 = delta_hat/(n*delta_hatnon)
        Var_varDLN2 = (sigmaSq_hat/size_nonzero) + (2*sigmaSq_hat^2/(size_nonzero-1))
        Var_varDLN = Var_varDLN1 + (4*Var_varDLN2)
        
        Var_varDLN12 = Var_varDLN[1]+Var_varDLN[2]
        Var_varDLN13 = Var_varDLN[1]+Var_varDLN[3]
        Var_varDLN14 = Var_varDLN[1]+Var_varDLN[4]
        Var_varDLN15 = Var_varDLN[1]+Var_varDLN[5]
        Var_varDLN23 = Var_varDLN[2]+Var_varDLN[3]
        Var_varDLN24 = Var_varDLN[2]+Var_varDLN[4]
        Var_varDLN25 = Var_varDLN[2]+Var_varDLN[5]
        Var_varDLN34 = Var_varDLN[3]+Var_varDLN[4]
        Var_varDLN35 = Var_varDLN[3]+Var_varDLN[5]
        Var_varDLN45 = Var_varDLN[4]+Var_varDLN[5]
        
        
        
        ##############################################################################
        #                 Algorithm 7: BR                                           #
        ##############################################################################   
        
        ############################### BR ######################################## 
        ## chain = 4, iter = 1250 to obtain m = 4*(1250/2) = 2500  
        ## refresh = 0 denoted as Stan running without output details
        x1 = as.numeric(dataframe.x1[rep,])
        x2 = as.numeric(dataframe.x2[rep,])
        x3 = as.numeric(dataframe.x3[rep,])
        x4 = as.numeric(dataframe.x4[rep,])
        x5 = as.numeric(dataframe.x5[rep,])
        
        fit.ref = stan('Stan_RBprior_K5.stan', data=list(n1=n[1], sizenonzero1=size_nonzero[1], 
                                                               y1=x1, sizezero1=size_zero[1], mu1=mu_hat[1], 
                                                               sigmaSq1=sigmaSq_hat[1], delta1=delta_hat[1],
                                                               
                                                               n2=n[2], sizenonzero2=size_nonzero[2], 
                                                               y2=x2, sizezero2=size_zero[2], mu2=mu_hat[2], 
                                                               sigmaSq2=sigmaSq_hat[2], delta2=delta_hat[2],
                                                               
                                                               n3=n[3], sizenonzero3=size_nonzero[3], 
                                                               y3=x3, sizezero3=size_zero[3], mu3=mu_hat[3], 
                                                               sigmaSq3=sigmaSq_hat[3], delta3=delta_hat[3],
                                                               
                                                               n4=n[4], sizenonzero4=size_nonzero[4], 
                                                               y4=x4, sizezero4=size_zero[4], mu4=mu_hat[4], 
                                                               sigmaSq4=sigmaSq_hat[4], delta4=delta_hat[4], 
                                                               
                                                               n5=n[5], sizenonzero5=size_nonzero[5], 
                                                               y5=x5, sizezero5=size_zero[5], mu5=mu_hat[5], 
                                                               sigmaSq5=sigmaSq_hat[5], delta5=delta_hat[5]),
                       iter = 1250, refresh = 0)
        
        ## Extracting the posterior draws
        post.mu.ref1 = extract(fit.ref)$mu1
        post.sigmaSq.ref1 = extract(fit.ref)$sigmaSq1
        post.delta.ref1 = extract(fit.ref)$delta1
        post_varDLN.ref1 = log(1-post.delta.ref1)+2*(post.mu.ref1+post.sigmaSq.ref1)
        
        post.mu.ref2 = extract(fit.ref)$mu2
        post.sigmaSq.ref2 = extract(fit.ref)$sigmaSq2
        post.delta.ref2 = extract(fit.ref)$delta2
        post_varDLN.ref2 = log(1-post.delta.ref2)+2*(post.mu.ref2+post.sigmaSq.ref2)
        
        post.mu.ref3 = extract(fit.ref)$mu3
        post.sigmaSq.ref3 = extract(fit.ref)$sigmaSq3
        post.delta.ref3 = extract(fit.ref)$delta3
        post_varDLN.ref3 = log(1-post.delta.ref3)+2*(post.mu.ref3+post.sigmaSq.ref3)
        
        post.mu.ref4 = extract(fit.ref)$mu4
        post.sigmaSq.ref4 = extract(fit.ref)$sigmaSq4
        post.delta.ref4 = extract(fit.ref)$delta4
        post_varDLN.ref4 = log(1-post.delta.ref4)+2*(post.mu.ref4+post.sigmaSq.ref4)
        
        post.mu.ref5 = extract(fit.ref)$mu5
        post.sigmaSq.ref5 = extract(fit.ref)$sigmaSq5
        post.delta.ref5 = extract(fit.ref)$delta5
        post_varDLN.ref5 = log(1-post.delta.ref5)+2*(post.mu.ref5+post.sigmaSq.ref5)
        
        post_varDLN.ref12 = post_varDLN.ref1-post_varDLN.ref2
        post_varDLN.ref13 = post_varDLN.ref1-post_varDLN.ref3
        post_varDLN.ref14 = post_varDLN.ref1-post_varDLN.ref4
        post_varDLN.ref15 = post_varDLN.ref1-post_varDLN.ref5
        post_varDLN.ref23 = post_varDLN.ref2-post_varDLN.ref3
        post_varDLN.ref24 = post_varDLN.ref2-post_varDLN.ref4
        post_varDLN.ref25 = post_varDLN.ref2-post_varDLN.ref5
        post_varDLN.ref34 = post_varDLN.ref3-post_varDLN.ref4
        post_varDLN.ref35 = post_varDLN.ref3-post_varDLN.ref5
        post_varDLN.ref45 = post_varDLN.ref4-post_varDLN.ref5
        
        for(idx.BR in 1:m){
          D12.BR = abs(thetahat.12 - post_varDLN.ref12[idx.BR])
          D13.BR = abs(thetahat.13 - post_varDLN.ref13[idx.BR])
          D14.BR = abs(thetahat.14 - post_varDLN.ref14[idx.BR])
          D15.BR = abs(thetahat.15 - post_varDLN.ref15[idx.BR])
          D23.BR = abs(thetahat.23 - post_varDLN.ref23[idx.BR])
          D24.BR = abs(thetahat.24 - post_varDLN.ref24[idx.BR])
          D25.BR = abs(thetahat.25 - post_varDLN.ref25[idx.BR])
          D34.BR = abs(thetahat.34 - post_varDLN.ref34[idx.BR])
          D35.BR = abs(thetahat.35 - post_varDLN.ref35[idx.BR])
          D45.BR = abs(thetahat.45 - post_varDLN.ref45[idx.BR])
          
          A.BR = max(D12.BR, D13.BR, D14.BR, D15.BR, D23.BR, D24.BR, D25.BR, D34.BR, D35.BR, D45.BR)
          B.BR = min(D12.BR, D13.BR, D14.BR, D15.BR, D23.BR, D24.BR, D25.BR, D34.BR, D35.BR, D45.BR)
          
          D.BR[idx.BR] = A.BR - B.BR
          
        }
        
        
        q.BR =   quantile(D.BR,0.95,type=8)
        L.BR12 = thetahat.12 - q.BR                     
        U.BR12 = thetahat.12 + q.BR
        CI.BR12 = data.frame(rep,"CI.BR12",L.BR12,U.BR12,theta12)
        names(CI.BR12) = c("rep","CI","lower","upper","theta")
        CI.BR = rbind(CI.BR,CI.BR12) 
        
        L.BR13 = thetahat.13 - q.BR                       
        U.BR13 = thetahat.13 + q.BR
        CI.BR13 = data.frame(rep,"CI.BR13",L.BR13,U.BR13,theta13)
        names(CI.BR13) = c("rep","CI","lower","upper","theta")
        CI.BR = rbind(CI.BR,CI.BR13) 
        
        L.BR14 = thetahat.14 - q.BR                       
        U.BR14 = thetahat.14 + q.BR
        CI.BR14 = data.frame(rep,"CI.BR14",L.BR14,U.BR14,theta14)
        names(CI.BR14) = c("rep","CI","lower","upper","theta")
        CI.BR = rbind(CI.BR,CI.BR14) 
        
        L.BR15 = thetahat.15 - q.BR                        
        U.BR15 = thetahat.15 + q.BR
        CI.BR15 = data.frame(rep,"CI.BR15",L.BR15,U.BR15,theta15)
        names(CI.BR15) = c("rep","CI","lower","upper","theta")
        CI.BR = rbind(CI.BR,CI.BR15) 
        
        L.BR23 = thetahat.23 - q.BR                    
        U.BR23 = thetahat.23 + q.BR
        CI.BR23 = data.frame(rep,"CI.BR23",L.BR23,U.BR23,theta23)
        names(CI.BR23) = c("rep","CI","lower","upper","theta")
        CI.BR = rbind(CI.BR,CI.BR23) 
        
        L.BR24 = thetahat.24 - q.BR                       
        U.BR24 = thetahat.24 + q.BR
        CI.BR24 = data.frame(rep,"CI.BR24",L.BR24,U.BR24,theta24)
        names(CI.BR24) = c("rep","CI","lower","upper","theta")
        CI.BR = rbind(CI.BR,CI.BR24) 
        
        L.BR25 = thetahat.25 - q.BR                        
        U.BR25 = thetahat.25 + q.BR
        CI.BR25 = data.frame(rep,"CI.BR25",L.BR25,U.BR25,theta25)
        names(CI.BR25) = c("rep","CI","lower","upper","theta")
        CI.BR = rbind(CI.BR,CI.BR25) 
        
        L.BR34 = thetahat.34 - q.BR                        
        U.BR34 = thetahat.34 + q.BR
        CI.BR34 = data.frame(rep,"CI.BR34",L.BR34,U.BR34,theta34)
        names(CI.BR34) = c("rep","CI","lower","upper","theta")
        CI.BR = rbind(CI.BR,CI.BR34) 
        
        L.BR35 = thetahat.35 - q.BR                        
        U.BR35 = thetahat.35 + q.BR
        CI.BR35 = data.frame(rep,"CI.BR35",L.BR35,U.BR35,theta35)
        names(CI.BR35) = c("rep","CI","lower","upper","theta")
        CI.BR = rbind(CI.BR,CI.BR35) 
        
        L.BR45 = thetahat.45 - q.BR                        
        U.BR45 = thetahat.45 + q.BR
        CI.BR45 = data.frame(rep,"CI.BR45",L.BR45,U.BR45,theta45)
        names(CI.BR45) = c("rep","CI","lower","upper","theta")
        CI.BR = rbind(CI.BR,CI.BR45)
        
        if(rep%%100 ==0){
          print(rep)
        }
      }
      
      CI.BR['CP.LE'] = 0
      CI.BR['CP'] = 0
      CI.BR['CP.UE'] = 0
      CI.BR[which(CI.BR["lower"]<CI.BR["theta"] & CI.BR["upper"]>CI.BR["theta"]),"CP"]=1
      CI.BR[which(CI.BR["lower"]>CI.BR["theta"]),"CP.LE"]=1
      CI.BR[which(CI.BR["upper"]<CI.BR["theta"]),"CP.UE"]=1
      CI.BR["Length"] = CI.BR["upper"] - CI.BR["lower"]
      ACI.BR = colMeans(CI.BR[,c("CP.LE","CP","CP.UE","Length")])
      ACP.BR = as.numeric(ACI.BR["CP"])*100
      ACP.LE.BR = as.numeric(ACI.BR["CP.LE"])*100
      ACP.UE.BR = as.numeric(ACI.BR["CP.UE"])*100
      Ex.length.BR = as.numeric(ACI.BR["Length"])	
      
      str1 = paste("M=",M,"m=",m, "mu=",paste(muTotal,collapse = ' '),"k=",k, "nTotal=",paste(nTotal,collapse = ' '))
      str2 = paste("sigmaSqTotal=", paste(sigmaSqTotal,collapse = ' '),"deltaTotal=",paste(deltaTotal,collapse = ' '),"theta=",paste(theta,collapse = ' '))   
      
      str3 =paste( "LE|BR =",paste(ACP.LE.BR,collapse = ' '),"CP|BR =", paste(ACP.BR,collapse = ' '),"UE|BR =", paste(ACP.UE.BR,collapse = ' '), "AL|BR =", paste(Ex.length.BR,collapse = ' '))
      
      alloutput = paste(str1,str2,str3,collapse = ' ')
      folderName = unlist(strsplit(currentFolder, "//"))
      dir.create(file.path(resultoutput), showWarnings = FALSE)
      write(alloutput, file = paste(resultoutput,"/Results_BR.txt",sep = ""), append = TRUE)
    }
    
    
    
    ##############################################################################
    #                     Algorithm 3: GPQ                                       #
    ##############################################################################
    resultoutput = "./GPQ_result/"
    
    #Delete all result
    delfiles = dir(path=resultoutput ,pattern="*.txt")
    file.remove(file.path(resultoutput, delfiles))
    
    M = 5000
    m = 2500
    
    
    listFolder = list.dirs("./DataSource/", recursive=FALSE)
    
    eachFolder = listFolder[1]
    print(listFolder)
    for (eachFolder in listFolder){
      
      print(paste("Start Folder",eachFolder,sep=""))
      ## read data from config file
      currentFolder = eachFolder
      configfile = read.csv(file=paste(currentFolder,"/config.csv",sep = ""))
      nTotal = configfile$config_ntotal
      sigmaSqTotal = configfile$config_sigmaSq
      deltaTotal = configfile$config_delta
      muTotal = -sigmaSqTotal/2
      k = length(nTotal)
      
      
      # calculate parameter from data in config file
      deltanonTotal = 1 - deltaTotal
      size_nonzeroTotal = nTotal*(1-deltaTotal)
      Ex_LN = exp(muTotal+(sigmaSqTotal/2))
      cv_LN = sqrt(exp(sigmaSqTotal)-1)
      theta = log(deltanonTotal)+(2*muTotal+sigmaSqTotal)+log(exp(sigmaSqTotal)-deltanonTotal)
      theta12 = theta[1]-theta[2]
      theta13 = theta[1]-theta[3]
      theta14 = theta[1]-theta[4]
      theta15 = theta[1]-theta[5]
      theta23 = theta[2]-theta[3]
      theta24 = theta[2]-theta[4]
      theta25 = theta[2]-theta[5]
      theta34 = theta[3]-theta[4]
      theta35 = theta[3]-theta[5]
      theta45 = theta[4]-theta[5]
      
      #read summary data
      summaryData = read.csv(file= paste(currentFolder,"/basicSummary_K",k,".csv",sep = ""))
      
      #GPQ
      CI.GPQ = data.frame()
      D.GPQ = rep(0,m)
      
      rep = 1
      for(rep in seq(1:M)){
        
        #Select data from repeat loop number
        currentData = summaryData[which(summaryData$loop==rep),]
        
        #Extract each parameter 
        delta_hatnon = currentData$delta_hatnon
        delta_hat = 1-delta_hatnon
        size_nonzero = currentData$size_non
        size_zero = currentData$size_zero
        mu_hat = currentData$mu_hat
        sigmaSq_hat = currentData$sigmaSq_hat
        n = size_nonzero + size_zero
        
        thetahat = log(delta_hatnon) + (2*(mu_hat+sigmaSq_hat))
        thetahat.12 = thetahat[1]-thetahat[2]
        thetahat.13 = thetahat[1]-thetahat[3]
        thetahat.14 = thetahat[1]-thetahat[4]
        thetahat.15 = thetahat[1]-thetahat[5]
        thetahat.23 = thetahat[2]-thetahat[3]
        thetahat.24 = thetahat[2]-thetahat[4]
        thetahat.25 = thetahat[2]-thetahat[5]
        thetahat.34 = thetahat[3]-thetahat[4]
        thetahat.35 = thetahat[3]-thetahat[5]
        thetahat.45 = thetahat[4]-thetahat[5]
        
        Var_varDLN1 = delta_hat/(n*delta_hatnon)
        Var_varDLN2 = (sigmaSq_hat/size_nonzero) + (2*sigmaSq_hat^2/(size_nonzero-1))
        Var_varDLN = Var_varDLN1 + (4*Var_varDLN2)
        
        Var_varDLN12 = Var_varDLN[1]+Var_varDLN[2]
        Var_varDLN13 = Var_varDLN[1]+Var_varDLN[3]
        Var_varDLN14 = Var_varDLN[1]+Var_varDLN[4]
        Var_varDLN15 = Var_varDLN[1]+Var_varDLN[5]
        Var_varDLN23 = Var_varDLN[2]+Var_varDLN[3]
        Var_varDLN24 = Var_varDLN[2]+Var_varDLN[4]
        Var_varDLN25 = Var_varDLN[2]+Var_varDLN[5]
        Var_varDLN34 = Var_varDLN[3]+Var_varDLN[4]
        Var_varDLN35 = Var_varDLN[3]+Var_varDLN[5]
        Var_varDLN45 = Var_varDLN[4]+Var_varDLN[5]
        
        
        
        #####################################################################
        #                     Algorithm 3: GPQ                             #
        #####################################################################
        idx.GPQ = 1
        
        while(idx.GPQ <= m){
          
          U2_1 = rchisq(1,size_nonzero[1]-1)/(size_nonzero[1]-1)
          U2_2 = rchisq(1,size_nonzero[2]-1)/(size_nonzero[2]-1)
          U2_3 = rchisq(1,size_nonzero[3]-1)/(size_nonzero[3]-1)
          U2_4 = rchisq(1,size_nonzero[4]-1)/(size_nonzero[4]-1)
          U2_5 = rchisq(1,size_nonzero[5]-1)/(size_nonzero[5]-1)
          U2 = c(U2_1, U2_2, U2_3, U2_4, U2_5)
          
          mu_GPQ1 = mu_hat[1] - ((rnorm(1)/sqrt(U2[1]))*sqrt(sigmaSq_hat[1]/size_nonzero[1]))
          mu_GPQ2 = mu_hat[2] - ((rnorm(1)/sqrt(U2[2]))*sqrt(sigmaSq_hat[2]/size_nonzero[2]))
          mu_GPQ3 = mu_hat[3] - ((rnorm(1)/sqrt(U2[3]))*sqrt(sigmaSq_hat[3]/size_nonzero[3]))
          mu_GPQ4 = mu_hat[4] - ((rnorm(1)/sqrt(U2[4]))*sqrt(sigmaSq_hat[4]/size_nonzero[4]))
          mu_GPQ5 = mu_hat[5] - ((rnorm(1)/sqrt(U2[5]))*sqrt(sigmaSq_hat[5]/size_nonzero[5]))
          mu_GPQ = c(mu_GPQ1, mu_GPQ2, mu_GPQ3, mu_GPQ4, mu_GPQ5)
          sigmaSq_GPQ = sigmaSq_hat/U2
          
          delta_GPQ1 = (sin(asin(sqrt(delta_hat[1]))-(rnorm(1)/(2*sqrt(n[1])))))^2
          delta_GPQ2 = (sin(asin(sqrt(delta_hat[2]))-(rnorm(1)/(2*sqrt(n[2])))))^2
          delta_GPQ3 = (sin(asin(sqrt(delta_hat[3]))-(rnorm(1)/(2*sqrt(n[3])))))^2
          delta_GPQ4 = (sin(asin(sqrt(delta_hat[4]))-(rnorm(1)/(2*sqrt(n[4])))))^2
          delta_GPQ5 = (sin(asin(sqrt(delta_hat[5]))-(rnorm(1)/(2*sqrt(n[5])))))^2
          delta_GPQ = c(delta_GPQ1, delta_GPQ2, delta_GPQ3, delta_GPQ4, delta_GPQ5)
          
          theta_GPQ_i  = log(1-delta_GPQ) + 2*(mu_GPQ+ sigmaSq_GPQ)
          theta_GPQ12 = theta_GPQ_i[1]-theta_GPQ_i[2]
          theta_GPQ13 = theta_GPQ_i[1]-theta_GPQ_i[3]
          theta_GPQ14 = theta_GPQ_i[1]-theta_GPQ_i[4]
          theta_GPQ15 = theta_GPQ_i[1]-theta_GPQ_i[5]
          theta_GPQ23 = theta_GPQ_i[2]-theta_GPQ_i[3]
          theta_GPQ24 = theta_GPQ_i[2]-theta_GPQ_i[4]
          theta_GPQ25 = theta_GPQ_i[2]-theta_GPQ_i[5]
          theta_GPQ34 = theta_GPQ_i[3]-theta_GPQ_i[4]
          theta_GPQ35 = theta_GPQ_i[3]-theta_GPQ_i[5]
          theta_GPQ45 = theta_GPQ_i[4]-theta_GPQ_i[5]
          
          D12.GPQ = abs((thetahat.12-theta_GPQ12)/sqrt(Var_varDLN12))
          D13.GPQ = abs((thetahat.13-theta_GPQ13)/sqrt(Var_varDLN13))
          D14.GPQ = abs((thetahat.14-theta_GPQ14)/sqrt(Var_varDLN14))
          D15.GPQ = abs((thetahat.15-theta_GPQ15)/sqrt(Var_varDLN15))
          D23.GPQ = abs((thetahat.23-theta_GPQ23)/sqrt(Var_varDLN23))
          D24.GPQ = abs((thetahat.24-theta_GPQ24)/sqrt(Var_varDLN24))
          D25.GPQ = abs((thetahat.25-theta_GPQ25)/sqrt(Var_varDLN25))
          D34.GPQ = abs((thetahat.34-theta_GPQ34)/sqrt(Var_varDLN34))
          D35.GPQ = abs((thetahat.35-theta_GPQ35)/sqrt(Var_varDLN35))
          D45.GPQ = abs((thetahat.45-theta_GPQ45)/sqrt(Var_varDLN45))
          
          D.GPQ[idx.GPQ] = max(D12.GPQ, D13.GPQ, D14.GPQ, D15.GPQ, D23.GPQ, D24.GPQ, D25.GPQ, D34.GPQ, D35.GPQ, D45.GPQ)
          idx.GPQ=idx.GPQ+1
        }
        
        q.GPQ =   quantile(D.GPQ,0.95,type=8)
        L.GPQ12 = thetahat.12-(q.GPQ*(sqrt(Var_varDLN12)))                         
        U.GPQ12 = thetahat.12+(q.GPQ*(sqrt(Var_varDLN12)))
        CI.GPQ12 = data.frame(rep,"CI.GPQ12",L.GPQ12,U.GPQ12,theta12)
        names(CI.GPQ12) = c("rep","CI","lower","upper","theta")
        CI.GPQ = rbind(CI.GPQ,CI.GPQ12)                                    
        
        L.GPQ13 = thetahat.13-(q.GPQ*(sqrt(Var_varDLN13)))                         
        U.GPQ13 = thetahat.13+(q.GPQ*(sqrt(Var_varDLN13)))
        CI.GPQ13 = data.frame(rep,"CI.GPQ13",L.GPQ13,U.GPQ13,theta13)
        names(CI.GPQ13) = c("rep","CI","lower","upper","theta")
        CI.GPQ = rbind(CI.GPQ,CI.GPQ13)                                        
        
        L.GPQ14 = thetahat.14-(q.GPQ*(sqrt(Var_varDLN14)))                         
        U.GPQ14 = thetahat.14+(q.GPQ*(sqrt(Var_varDLN14)))
        CI.GPQ14 = data.frame(rep,"CI.GPQ14",L.GPQ14,U.GPQ14,theta14)
        names(CI.GPQ14) = c("rep","CI","lower","upper","theta")
        CI.GPQ = rbind(CI.GPQ,CI.GPQ14)                                          
        
        L.GPQ15 = thetahat.15-(q.GPQ*(sqrt(Var_varDLN15)))                         
        U.GPQ15 = thetahat.15+(q.GPQ*(sqrt(Var_varDLN15)))
        CI.GPQ15 = data.frame(rep,"CI.GPQ15",L.GPQ15,U.GPQ15,theta15)
        names(CI.GPQ15) = c("rep","CI","lower","upper","theta")
        CI.GPQ = rbind(CI.GPQ,CI.GPQ15)                                  
        
        L.GPQ23 = thetahat.23-(q.GPQ*(sqrt(Var_varDLN23)))                         
        U.GPQ23 = thetahat.23+(q.GPQ*(sqrt(Var_varDLN23)))
        CI.GPQ23 = data.frame(rep,"CI.GPQ23",L.GPQ23,U.GPQ23,theta23)
        names(CI.GPQ23) = c("rep","CI","lower","upper","theta")
        CI.GPQ = rbind(CI.GPQ,CI.GPQ23)                                       
        
        L.GPQ24 = thetahat.24-(q.GPQ*(sqrt(Var_varDLN24)))                        
        U.GPQ24 = thetahat.24+(q.GPQ*(sqrt(Var_varDLN24)))
        CI.GPQ24 = data.frame(rep,"CI.GPQ24",L.GPQ24,U.GPQ24,theta24)
        names(CI.GPQ24) = c("rep","CI","lower","upper","theta")
        CI.GPQ = rbind(CI.GPQ,CI.GPQ24)                                       
        
        L.GPQ25 = thetahat.25-(q.GPQ*(sqrt(Var_varDLN25)))                         
        U.GPQ25 = thetahat.25+(q.GPQ*(sqrt(Var_varDLN25)))
        CI.GPQ25 = data.frame(rep,"CI.GPQ25",L.GPQ25,U.GPQ25,theta25)
        names(CI.GPQ25) = c("rep","CI","lower","upper","theta")
        CI.GPQ = rbind(CI.GPQ,CI.GPQ25)                                        
        
        L.GPQ34 = thetahat.34-(q.GPQ*(sqrt(Var_varDLN34)))                         
        U.GPQ34 = thetahat.34+(q.GPQ*(sqrt(Var_varDLN34)))
        CI.GPQ34 = data.frame(rep,"CI.GPQ34",L.GPQ34,U.GPQ34,theta34)
        names(CI.GPQ34) = c("rep","CI","lower","upper","theta")
        CI.GPQ = rbind(CI.GPQ,CI.GPQ34)                                        
        
        L.GPQ35 = thetahat.35-(q.GPQ*(sqrt(Var_varDLN35)))                         
        U.GPQ35 = thetahat.35+(q.GPQ*(sqrt(Var_varDLN35)))
        CI.GPQ35 = data.frame(rep,"CI.GPQ35",L.GPQ35,U.GPQ35,theta35)
        names(CI.GPQ35) = c("rep","CI","lower","upper","theta")
        CI.GPQ = rbind(CI.GPQ,CI.GPQ35)                                        
        
        L.GPQ45 = thetahat.45-(q.GPQ*(sqrt(Var_varDLN45)))                         
        U.GPQ45 = thetahat.45+(q.GPQ*(sqrt(Var_varDLN45)))
        CI.GPQ45 = data.frame(rep,"CI.GPQ45",L.GPQ45,U.GPQ45,theta45)
        names(CI.GPQ45) = c("rep","CI","lower","upper","theta")
        CI.GPQ = rbind(CI.GPQ,CI.GPQ45)  
        
      }
      
      
      CI.GPQ['CP.LE'] = 0
      CI.GPQ['CP'] = 0
      CI.GPQ['CP.UE'] = 0
      CI.GPQ[which(CI.GPQ["lower"]<CI.GPQ["theta"] & CI.GPQ["upper"]>CI.GPQ["theta"]),"CP"]=1
      CI.GPQ[which(CI.GPQ["lower"]>CI.GPQ["theta"]),"CP.LE"]=1
      CI.GPQ[which(CI.GPQ["upper"]<CI.GPQ["theta"]),"CP.UE"]=1
      CI.GPQ["Length"] = CI.GPQ["upper"] - CI.GPQ["lower"]
      ACI.GPQ = colMeans(CI.GPQ[,c("CP.LE","CP","CP.UE","Length")])
      ACP.GPQ = as.numeric(ACI.GPQ["CP"])*100
      ACP.LE.GPQ = as.numeric(ACI.GPQ["CP.LE"])*100
      ACP.UE.GPQ = as.numeric(ACI.GPQ["CP.UE"])*100
      Ex.length.GPQ = as.numeric(ACI.GPQ["Length"])	  
      
      str1 = paste("M=",M,"m=",m, "mu=",paste(muTotal,collapse = ' '),"k=",k, "nTotal=",paste(nTotal,collapse = ' '))
      str2 = paste("sigmaSqTotal=", paste(sigmaSqTotal,collapse = ' '),"deltaTotal=",paste(deltaTotal,collapse = ' '),"theta=",paste(theta,collapse = ' '))   
      
      str3 =paste( "LE|GPQ =",paste(ACP.LE.GPQ,collapse = ' '),"CP|GPQ =", paste(ACP.GPQ,collapse = ' '),"UE|GPQ =", paste(ACP.UE.GPQ,collapse = ' '), "AL|GPQ =", paste(Ex.length.GPQ,collapse = ' '))
      
      alloutput = paste(str1,str2,str3,collapse = ' ')
      folderName = unlist(strsplit(currentFolder, "//"))
      dir.create(file.path(resultoutput), showWarnings = FALSE)
      write(alloutput, file = paste(resultoutput,"/Results_GPQ.txt",sep = ""), append = TRUE)
    }
    
    
    
    ##############################################################################
    #                  Algorithm 4: PB                                           #
    ##############################################################################   
    ### This function for loading the raw dataset
    loadXlist = function(folderselect,Kselect){
      D1 = read.csv(file=paste(folderselect,"/D",Kselect,"_xlist.csv",sep = ""),header = FALSE)
      return (D1)
    }
    
    
    resultoutput = "./PB_result/"
    
    #Delete all result
    delfiles = dir(path=resultoutput ,pattern="*.txt")
    file.remove(file.path(resultoutput, delfiles))
    
    M = 5000
    m = 2500
    
    
    listFolder = list.dirs("./DataSource/", recursive=FALSE)
    
    eachFolder = listFolder[1]
    print(listFolder)
    for (eachFolder in listFolder){
      
      print(paste("Start Folder",eachFolder,sep=""))
      ## read data from config file
      currentFolder = eachFolder
      configfile = read.csv(file=paste(currentFolder,"/config.csv",sep = ""))
      nTotal = configfile$config_ntotal
      sigmaSqTotal = configfile$config_sigmaSq
      deltaTotal = configfile$config_delta
      muTotal = -sigmaSqTotal/2
      k = length(nTotal)
      
      
      # calculate parameter from data in config file
      deltanonTotal = 1 - deltaTotal
      size_nonzeroTotal = nTotal*(1-deltaTotal)
      Ex_LN = exp(muTotal+(sigmaSqTotal/2))
      cv_LN = sqrt(exp(sigmaSqTotal)-1)
      theta = log(deltanonTotal)+(2*muTotal+sigmaSqTotal)+log(exp(sigmaSqTotal)-deltanonTotal)
      theta12 = theta[1]-theta[2]
      theta13 = theta[1]-theta[3]
      theta14 = theta[1]-theta[4]
      theta15 = theta[1]-theta[5]
      theta23 = theta[2]-theta[3]
      theta24 = theta[2]-theta[4]
      theta25 = theta[2]-theta[5]
      theta34 = theta[3]-theta[4]
      theta35 = theta[3]-theta[5]
      theta45 = theta[4]-theta[5]
      
      #read summary data
      summaryData = read.csv(file= paste(currentFolder,"/basicSummary_K",k,".csv",sep = ""))
      
      #PB: Parametric bootstrap
      CI.PB = data.frame()
      D.PB = rep(0,m)
      
      rep = 1
      dataframe.x1 = loadXlist(currentFolder,1)
      dataframe.x2 = loadXlist(currentFolder,2)
      dataframe.x3 = loadXlist(currentFolder,3)
      dataframe.x4 = loadXlist(currentFolder,4)
      dataframe.x5 = loadXlist(currentFolder,5)
      for(rep in seq(1:M)){
        
        #Select data from repeat loop number
        currentData = summaryData[which(summaryData$loop==rep),]
        
        #Extract each parameter 
        delta_hatnon = currentData$delta_hatnon
        delta_hat = 1-delta_hatnon
        size_nonzero = currentData$size_non
        size_zero = currentData$size_zero
        mu_hat = currentData$mu_hat
        sigmaSq_hat = currentData$sigmaSq_hat
        n = size_nonzero + size_zero
        
        thetahat = log(delta_hatnon) + (2*(mu_hat+sigmaSq_hat))
        thetahat.12 = thetahat[1]-thetahat[2]
        thetahat.13 = thetahat[1]-thetahat[3]
        thetahat.14 = thetahat[1]-thetahat[4]
        thetahat.15 = thetahat[1]-thetahat[5]
        thetahat.23 = thetahat[2]-thetahat[3]
        thetahat.24 = thetahat[2]-thetahat[4]
        thetahat.25 = thetahat[2]-thetahat[5]
        thetahat.34 = thetahat[3]-thetahat[4]
        thetahat.35 = thetahat[3]-thetahat[5]
        thetahat.45 = thetahat[4]-thetahat[5]
        
        Var_varDLN1 = delta_hat/(n*delta_hatnon)
        Var_varDLN2 = (sigmaSq_hat/size_nonzero) + (2*sigmaSq_hat^2/(size_nonzero-1))
        Var_varDLN = Var_varDLN1 + (4*Var_varDLN2)
        
        Var_varDLN12 = Var_varDLN[1]+Var_varDLN[2]
        Var_varDLN13 = Var_varDLN[1]+Var_varDLN[3]
        Var_varDLN14 = Var_varDLN[1]+Var_varDLN[4]
        Var_varDLN15 = Var_varDLN[1]+Var_varDLN[5]
        Var_varDLN23 = Var_varDLN[2]+Var_varDLN[3]
        Var_varDLN24 = Var_varDLN[2]+Var_varDLN[4]
        Var_varDLN25 = Var_varDLN[2]+Var_varDLN[5]
        Var_varDLN34 = Var_varDLN[3]+Var_varDLN[4]
        Var_varDLN35 = Var_varDLN[3]+Var_varDLN[5]
        Var_varDLN45 = Var_varDLN[4]+Var_varDLN[5]
        
        
        
        ##############################################################################
        #                  Algorithm 4: PB interval                                  #
        ##############################################################################    
        x1 = as.numeric(dataframe.x1[rep,])
        x2 = as.numeric(dataframe.x2[rep,])
        x3 = as.numeric(dataframe.x3[rep,])
        x4 = as.numeric(dataframe.x4[rep,])
        x5 = as.numeric(dataframe.x5[rep,])
        
        idx.PB = 1
        
        while(idx.PB <= m){
          
          #Bootstrap
          A1 = sample(1:n[1], size=n[1], replace=TRUE)
          X1.BS = x1[A1]
          nonzero_value1.BS = log(X1.BS[which(X1.BS!=0)])                          #find values is not zero
          size_nonzero1.BS =length(nonzero_value1.BS) 
          size_zero1.BS = n[1]-size_nonzero1.BS
          mu_hat1.BS = mean(nonzero_value1.BS)
          sigmaSq_hat1.BS = var(nonzero_value1.BS)
          delta_hat1.BS = 1-(size_nonzero1.BS/n[1])  
          
          A2 = sample(1:n[2], size=n[2], replace=TRUE)
          X2.BS = x2[A2]
          nonzero_value2.BS = log(X2.BS[which(X2.BS!=0)])                          #find values is not zero
          size_nonzero2.BS =length(nonzero_value2.BS) 
          size_zero2.BS = n[2]-size_nonzero2.BS
          mu_hat2.BS = mean(nonzero_value2.BS)
          sigmaSq_hat2.BS = var(nonzero_value2.BS)
          delta_hat2.BS = 1-(size_nonzero2.BS/n[2])  
          
          A3 = sample(1:n[3], size=n[3], replace=TRUE)
          X3.BS = x3[A3]
          nonzero_value3.BS = log(X3.BS[which(X3.BS!=0)])                          #find values is not zero
          size_nonzero3.BS =length(nonzero_value3.BS) 
          size_zero3.BS = n[3]-size_nonzero3.BS
          mu_hat3.BS = mean(nonzero_value3.BS)
          sigmaSq_hat3.BS = var(nonzero_value3.BS)
          delta_hat3.BS = 1-(size_nonzero3.BS/n[3])
          
          A4 = sample(1:n[4], size=n[4], replace=TRUE)
          X4.BS = x4[A4]
          nonzero_value4.BS = log(X4.BS[which(X4.BS!=0)])                          #find values is not zero
          size_nonzero4.BS =length(nonzero_value4.BS) 
          size_zero4.BS = n[4]-size_nonzero4.BS
          mu_hat4.BS = mean(nonzero_value4.BS)
          sigmaSq_hat4.BS = var(nonzero_value4.BS)
          delta_hat4.BS = 1-(size_nonzero4.BS/n[4]) 
          
          A5 =sample(1:n[5],size=n[5],replace=TRUE)
          X5.BS = x5[A5]
          nonzero_value5.BS = log(X5.BS[which(X5.BS!=0)])                          
          size_nonzero5.BS =length(nonzero_value5.BS) 
          size_zero5.BS = n[5]-size_nonzero5.BS
          mu_hat5.BS = mean(nonzero_value5.BS)
          sigmaSq_hat5.BS = var(nonzero_value5.BS)
          delta_hat5.BS = 1-(size_nonzero5.BS/n[5])
          
          
          size_nonzero.BS = c(size_nonzero1.BS, size_nonzero2.BS, size_nonzero3.BS, size_nonzero4.BS, size_nonzero5.BS)
          mu_hat.BS = c(mu_hat1.BS, mu_hat2.BS, mu_hat3.BS, mu_hat4.BS, mu_hat5.BS)
          sigmaSq_hat.BS = c(sigmaSq_hat1.BS, sigmaSq_hat2.BS, sigmaSq_hat3.BS, sigmaSq_hat4.BS, sigmaSq_hat5.BS)
          delta_hat.BS = c(delta_hat1.BS, delta_hat2.BS, delta_hat3.BS, delta_hat4.BS, delta_hat5.BS)
          delta_hatnon.BS = 1 - delta_hat.BS
          thetahat.BS = log(delta_hatnon.BS) + (2*(mu_hat.BS+sigmaSq_hat.BS))
          
          
          #Parametric boostrap
          mu_hat1.PB =  mu_hat1.BS + (rnorm(1)*sqrt(sigmaSq_hat1.BS/size_nonzero1.BS))
          mu_hat2.PB =  mu_hat2.BS + (rnorm(1)*sqrt(sigmaSq_hat2.BS/size_nonzero2.BS))
          mu_hat3.PB =  mu_hat3.BS + (rnorm(1)*sqrt(sigmaSq_hat3.BS/size_nonzero3.BS))
          mu_hat4.PB =  mu_hat4.BS + (rnorm(1)*sqrt(sigmaSq_hat4.BS/size_nonzero4.BS))
          mu_hat5.PB =  mu_hat5.BS + (rnorm(1)*sqrt(sigmaSq_hat5.BS/size_nonzero5.BS))
          
          sigmaSq_hat1.PB = sigmaSq_hat1.BS*rchisq(1,size_nonzero1.BS-1, ncp = 0)/(size_nonzero1.BS-1)
          sigmaSq_hat2.PB = sigmaSq_hat2.BS*rchisq(1,size_nonzero2.BS-1, ncp = 0)/(size_nonzero2.BS-1)
          sigmaSq_hat3.PB = sigmaSq_hat3.BS*rchisq(1,size_nonzero3.BS-1, ncp = 0)/(size_nonzero3.BS-1)
          sigmaSq_hat4.PB = sigmaSq_hat4.BS*rchisq(1,size_nonzero4.BS-1, ncp = 0)/(size_nonzero4.BS-1)
          sigmaSq_hat5.PB = sigmaSq_hat5.BS*rchisq(1,size_nonzero5.BS-1, ncp = 0)/(size_nonzero5.BS-1)
          
          delta_hat1.PB = rbeta(1, size_zero1.BS+0.5, size_nonzero1.BS+0.5)
          delta_hat2.PB = rbeta(1, size_zero2.BS+0.5, size_nonzero2.BS+0.5)
          delta_hat3.PB = rbeta(1, size_zero3.BS+0.5, size_nonzero3.BS+0.5)
          delta_hat4.PB = rbeta(1, size_zero4.BS+0.5, size_nonzero4.BS+0.5)
          delta_hat5.PB = rbeta(1, size_zero5.BS+0.5, size_nonzero5.BS+0.5)
          delta_hat.PB = c(delta_hat1.PB, delta_hat2.PB, delta_hat3.PB, delta_hat4.PB, delta_hat5.PB)
          delta_hatnon.PB = 1-delta_hat.PB
          size_nonzero.PB = n*delta_hatnon.PB
          mu_hat.PB = c(mu_hat1.PB, mu_hat2.PB, mu_hat3.PB, mu_hat4.PB, mu_hat5.PB)
          sigmaSq_hat.PB = c(sigmaSq_hat1.PB, sigmaSq_hat2.PB, sigmaSq_hat3.PB, sigmaSq_hat4.PB, sigmaSq_hat5.PB)
          thetahat.PB = log(delta_hatnon.PB) + (2*(mu_hat.PB+sigmaSq_hat.PB))
          
          Var_varDLN1.PB = delta_hat.PB/(n*delta_hatnon.PB)
          Var_varDLN2.PB = (sigmaSq_hat.PB/size_nonzero.PB) + (2*sigmaSq_hat.PB^2/(size_nonzero.PB-1))
          Var_varDLN.PB = Var_varDLN1.PB + (4*Var_varDLN2.PB) 
          V12.PB = Var_varDLN.PB[1] + Var_varDLN.PB[2]
          V13.PB = Var_varDLN.PB[1] + Var_varDLN.PB[3]
          V14.PB = Var_varDLN.PB[1] + Var_varDLN.PB[4]
          V15.PB = Var_varDLN.PB[1] + Var_varDLN.PB[5]
          V23.PB = Var_varDLN.PB[2] + Var_varDLN.PB[3]
          V24.PB = Var_varDLN.PB[2] + Var_varDLN.PB[4]
          V25.PB = Var_varDLN.PB[2] + Var_varDLN.PB[5]
          V34.PB = Var_varDLN.PB[3] + Var_varDLN.PB[4]
          V35.PB = Var_varDLN.PB[3] + Var_varDLN.PB[5]
          V45.PB = Var_varDLN.PB[4] + Var_varDLN.PB[5]
          
          
          D12.PB  = abs(((thetahat.PB[1] - thetahat.PB[2])-(thetahat.BS[1] - thetahat.BS[2]))/sqrt(V12.PB))
          D13.PB  = abs(((thetahat.PB[1] - thetahat.PB[3])-(thetahat.BS[1] - thetahat.BS[3]))/sqrt(V13.PB))
          D14.PB  = abs(((thetahat.PB[1] - thetahat.PB[4])-(thetahat.BS[1] - thetahat.BS[4]))/sqrt(V14.PB))
          D15.PB  = abs(((thetahat.PB[1] - thetahat.PB[5])-(thetahat.BS[1] - thetahat.BS[5]))/sqrt(V15.PB))
          D23.PB  = abs(((thetahat.PB[2] - thetahat.PB[3])-(thetahat.BS[2] - thetahat.BS[3]))/sqrt(V23.PB))
          D24.PB  = abs(((thetahat.PB[2] - thetahat.PB[4])-(thetahat.BS[2] - thetahat.BS[4]))/sqrt(V24.PB))
          D25.PB  = abs(((thetahat.PB[2] - thetahat.PB[5])-(thetahat.BS[2] - thetahat.BS[5]))/sqrt(V25.PB))
          D34.PB  = abs(((thetahat.PB[3] - thetahat.PB[4])-(thetahat.BS[3] - thetahat.BS[4]))/sqrt(V34.PB))
          D35.PB  = abs(((thetahat.PB[3] - thetahat.PB[5])-(thetahat.BS[3] - thetahat.BS[5]))/sqrt(V35.PB))
          D45.PB  = abs(((thetahat.PB[4] - thetahat.PB[5])-(thetahat.BS[4] - thetahat.BS[5]))/sqrt(V45.PB))
          
          D.PB[idx.PB] = max(D12.PB,D13.PB,D14.PB,D15.PB,D23.PB,D24.PB,D25.PB,D34.PB,D35.PB,D45.PB)
          if(is.na(D.PB[idx.PB])){
            next
          }
          
          idx.PB=idx.PB+1
          
        }
        
        q.PB = quantile(D.PB,0.95,type=8)
        L.PB12 = thetahat.12-(q.PB*(sqrt(Var_varDLN12)))                        
        U.PB12 = thetahat.12+(q.PB*(sqrt(Var_varDLN12)))  
        CI.PB12 = data.frame(rep,"CI.PB12",L.PB12,U.PB12,theta12)
        names(CI.PB12) = c("rep","CI","lower","upper","theta")
        CI.PB = rbind(CI.PB,CI.PB12)
        
        L.PB13 = thetahat.13-(q.PB*(sqrt(Var_varDLN13)))                         
        U.PB13 = thetahat.13+(q.PB*(sqrt(Var_varDLN13)))                        
        CI.PB13 = data.frame(rep,"CI.PB13",L.PB13,U.PB13,theta13)
        names(CI.PB13) = c("rep","CI","lower","upper","theta")
        CI.PB = rbind(CI.PB,CI.PB13)                                   
        
        L.PB14 = thetahat.14-(q.PB*(sqrt(Var_varDLN14)))                         
        U.PB14 = thetahat.14+(q.PB*(sqrt(Var_varDLN14)))                         
        CI.PB14 = data.frame(rep,"CI.PB14",L.PB14,U.PB14,theta14)
        names(CI.PB14) = c("rep","CI","lower","upper","theta")
        CI.PB = rbind(CI.PB,CI.PB14)                                     
        
        L.PB15 = thetahat.15-(q.PB*(sqrt(Var_varDLN15)))                        
        U.PB15 = thetahat.15+(q.PB*(sqrt(Var_varDLN15)))                         
        CI.PB15 = data.frame(rep,"CI.PB15",L.PB15,U.PB15,theta15)
        names(CI.PB15) = c("rep","CI","lower","upper","theta")
        CI.PB = rbind(CI.PB,CI.PB15)                                   
        
        L.PB23 = thetahat.23-(q.PB*(sqrt(Var_varDLN23)))                         
        U.PB23 = thetahat.23+(q.PB*(sqrt(Var_varDLN23)))                         
        CI.PB23 = data.frame(rep,"CI.PB23",L.PB23,U.PB23,theta23)
        names(CI.PB23) = c("rep","CI","lower","upper","theta")
        CI.PB = rbind(CI.PB,CI.PB23)                                    
        
        L.PB24 = thetahat.24-(q.PB*(sqrt(Var_varDLN24)))                         
        U.PB24 = thetahat.24+(q.PB*(sqrt(Var_varDLN24)))                         
        CI.PB24 = data.frame(rep,"CI.PB24",L.PB24,U.PB24,theta24)
        names(CI.PB24) = c("rep","CI","lower","upper","theta")
        CI.PB = rbind(CI.PB,CI.PB24)                                    
        
        L.PB25 = thetahat.25-(q.PB*(sqrt(Var_varDLN25)))                         
        U.PB25 = thetahat.25+(q.PB*(sqrt(Var_varDLN25)))                        
        CI.PB25 = data.frame(rep,"CI.PB25",L.PB25,U.PB25,theta25)
        names(CI.PB25) = c("rep","CI","lower","upper","theta")
        CI.PB = rbind(CI.PB,CI.PB25)                                    
        
        L.PB34 = thetahat.34-(q.PB*(sqrt(Var_varDLN34)))                        
        U.PB34 = thetahat.34+(q.PB*(sqrt(Var_varDLN34)))                         
        CI.PB34 = data.frame(rep,"CI.PB34",L.PB34,U.PB34,theta34)
        names(CI.PB34) = c("rep","CI","lower","upper","theta")
        CI.PB = rbind(CI.PB,CI.PB34)                                     
        
        L.PB35 = thetahat.35-(q.PB*(sqrt(Var_varDLN35)))                         
        U.PB35 = thetahat.35+(q.PB*(sqrt(Var_varDLN35)))                         
        CI.PB35 = data.frame(rep,"CI.PB35",L.PB35,U.PB35,theta35)
        names(CI.PB35) = c("rep","CI","lower","upper","theta")
        CI.PB = rbind(CI.PB,CI.PB35)                                     
        
        L.PB45 = thetahat.45-(q.PB*(sqrt(Var_varDLN45)))                        
        U.PB45 = thetahat.45+(q.PB*(sqrt(Var_varDLN45)))                         
        CI.PB45 = data.frame(rep,"CI.PB45",L.PB45,U.PB45,theta45)
        names(CI.PB45) = c("rep","CI","lower","upper","theta")
        CI.PB = rbind(CI.PB,CI.PB45)
        
        
      }
      
      CI.PB['CP.LE'] = 0
      CI.PB['CP'] = 0
      CI.PB['CP.UE'] = 0
      CI.PB[which(CI.PB["lower"]<CI.PB["theta"] & CI.PB["upper"]>CI.PB["theta"]),"CP"]=1
      CI.PB[which(CI.PB["lower"]>CI.PB["theta"]),"CP.LE"]=1
      CI.PB[which(CI.PB["upper"]<CI.PB["theta"]),"CP.UE"]=1
      CI.PB["Length"] = CI.PB["upper"] - CI.PB["lower"]
      ACI.PB = colMeans(CI.PB[,c("CP.LE","CP","CP.UE","Length")])
      ACP.PB = as.numeric(ACI.PB["CP"])*100
      ACP.LE.PB = as.numeric(ACI.PB["CP.LE"])*100
      ACP.UE.PB = as.numeric(ACI.PB["CP.UE"])*100
      Ex.length.PB = as.numeric(ACI.PB["Length"])	
      
      str1 = paste("M=",M,"m=",m, "mu=",paste(muTotal,collapse = ' '),"k=",k, "nTotal=",paste(nTotal,collapse = ' '))
      str2 = paste("sigmaSqTotal=", paste(sigmaSqTotal,collapse = ' '),"deltaTotal=",paste(deltaTotal,collapse = ' '),"theta=",paste(theta,collapse = ' '))   
      
      str3 =paste( "LE|PB =",paste(ACP.LE.PB,collapse = ' '),"CP|PB =", paste(ACP.PB,collapse = ' '),"UE|PB =", paste(ACP.UE.PB,collapse = ' '), "AL|PB =", paste(Ex.length.PB,collapse = ' '))
      
      alloutput = paste(str1,str2,str3,collapse = ' ')
      folderName = unlist(strsplit(currentFolder, "//"))
      dir.create(file.path(resultoutput), showWarnings = FALSE)
      write(alloutput, file = paste(resultoutput,"/Results_PB.txt",sep = ""), append = TRUE)
    }
    
    
    
    
    ##############################################################################
    # Daily precipitation Datasets                                               #
    ##############################################################################
    rm(list=ls())
    library(readxl)
    library(FSA)
    library(MASS)         #For fitdistr
    library(fitdistrplus)
    library(vcd)          #for expo distribution
    library(gridExtra)    #for multiplots
    
    #Dataset 1 : Songkhla station : SK
    SK = c(160.0,  14.6,  20.8,   8.8,   0.0,   0.0,   0.0,   0.2,  52.0,  39.4,   0.6,  12.2,   5.4,   9.4,   7.0,  19.2,  84.4,  97.2,  92.0,  19.8,
           5.4,   0.0,  23.8,  23.4,   2.2,   1.0,   0.0,   0.0,   0.0,   0.0,   6.2,   0.8,   1.4,   2.6,  21.4,   9.2,   0.2,   0.0,   0.4,   0.8,
           29.0,  23.0,   5.0,   0.0,   5.4,   1.8,   0.8,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,
           4.4,   9.6)
    
    
    #Dataset 2 : Songkhla Sadao station : SKSD
    SKSD = c(56.4, 85.8,  4.2,  0.2,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  4.2, 37.2,  0.0,  2.4, 25.6, 97.4,  9.2, 19.2,  7.2,  0.4,  0.0,  0.0,  0.0,  0.0, 10.0,  0.0,  0.0,
             0.0,  0.0,  0.4,  4.2,  8.2, 42.6,  8.4, 70.2,  2.8,  0.0,  0.0,  0.0, 15.6,  0.6,  0.2,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,
             0.0,  0.0,  0.0,  0.0,  0.0,  2.0)
    
     
    #Dataset 3 : Yala station : YL
    YL = c(46.4,  46.6,  55.8,  27.0,   0.2,   0.0,   0.0,   1.6,   0.0,   0.0,   2.8,  17.2,   2.0,   0.0,  12.4,  43.8, 126.4, 113.8,  39.8,  27.8,   0.0,   1.2,  31.0,
           19.6,  46.6,  27.6,   1.0,   0.0,   0.0,   0.0,  11.2,   6.6,   5.6,  49.6,  28.6, 137.8,  84.8,   1.8,   0.4,   0.0,   2.8,   0.2,   0.6,   2.4,   0.0,   0.0,
           0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   2.2,   0.0,   0.0,   2.0,   0.0,   0.0,   0.0,   3.0)
    
    
    #Dataset 4 : Narathiwat station : NRTW
    NRTW = c(38.6,  70.0,  74.2,   0.4,   0.0,   0.0,   0.0,   0.0,   0.0,   0.8,   9.2,   0.0,   8.2,   0.0,  78.4,  43.0, 162.0, 141.2,  43.6,  20.4,   0.2,   1.0,  61.4,
             6.6,  39.8,  84.0,   0.0,   0.0,   0.0,   0.0,  89.2,  31.2,   6.4,  38.6,  10.4,  62.8,  13.2,   9.2,   0.0,   1.4,  12.6,   0.2,   3.6,   3.0,   0.0,   0.0,
             0.0,   0.0,   1.2,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.2,   0.0,   0.0,   0.0,   0.4)
    
    
    #Dataset 5 : Pattani station : PTN
    PTN = c(82.0,   0.0,   0.0,   7.2,   0.0,   0.2,   0.0,   0.0,   0.0,   3.6,   9.8,   8.0,  12.8,   3.4,   7.2,  62.8, 164.8,  46.4,  26.2,   7.0,   3.4,   3.4,  12.6,
            0.0,   6.8,   2.8,   0.2,   0.0,   0.0,   0.0,   3.2,   0.8,   2.0,  49.8,   4.4,  49.0,   0.2,   2.8,   1.2,   0.0,  22.8,   0.0,   1.2,   1.0,   0.0,   0.0,
            0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   1.6)


    
    ##############################################################################
    # Basic statistics                                                            #
    ##############################################################################
    
    n_SK = length(SK)
    nonZero_SK = SK[which(SK!=0)]
    nonzero_value_SK = log(SK[which(SK!=0)])
    size_nonzero_SK = length(nonZero_SK)
    size_zero_SK = n_SK-length(nonZero_SK)
    mu_hat_SK = mean(nonzero_value_SK)
    sigmaSq_hat_SK  = var(nonzero_value_SK)
    delta_hatnon_SK = size_nonzero_SK/n_SK
    delta_hat_SK  = 1-delta_hatnon_SK 
    theta_hat_SK = log(delta_hatnon_SK) + 2*(mu_hat_SK+sigmaSq_hat_SK)
    

    n_SKSD = length(SKSD)
    nonZero_SKSD = SKSD[which(SKSD!=0)]
    nonzero_value_SKSD = log(SKSD[which(SKSD!=0)])
    size_nonzero_SKSD = length(nonZero_SKSD)
    size_zero_SKSD = n_SKSD-length(nonZero_SKSD)
    mu_hat_SKSD = mean(nonzero_value_SKSD)
    sigmaSq_hat_SKSD  = var(nonzero_value_SKSD)
    delta_hatnon_SKSD = size_nonzero_SKSD/n_SKSD
    delta_hat_SKSD  = 1-delta_hatnon_SKSD 
    theta_hat_SKSD = log(delta_hatnon_SKSD) + 2*(mu_hat_SKSD+sigmaSq_hat_SKSD)
    
    n_YL = length(YL)
    nonZero_YL = YL[which(YL!=0)]
    nonzero_value_YL = log(YL[which(YL!=0)])
    size_nonzero_YL = length(nonZero_YL)
    size_zero_YL = n_YL-length(nonZero_YL)
    mu_hat_YL = mean(nonzero_value_YL)
    sigmaSq_hat_YL  = var(nonzero_value_YL)
    delta_hatnon_YL = size_nonzero_YL/n_YL
    delta_hat_YL  = 1-delta_hatnon_YL 
    theta_hat_YL = log(delta_hatnon_YL) + 2*(mu_hat_YL+sigmaSq_hat_YL)
    

    n_NRTW = length(NRTW)
    nonZero_NRTW = NRTW[which(NRTW!=0)]
    nonzero_value_NRTW = log(NRTW[which(NRTW!=0)])
    size_nonzero_NRTW = length(nonZero_NRTW)
    size_zero_NRTW = n_NRTW-length(nonZero_NRTW)
    mu_hat_NRTW = mean(nonzero_value_NRTW)
    sigmaSq_hat_NRTW  = var(nonzero_value_NRTW)
    delta_hatnon_NRTW = size_nonzero_NRTW/n_NRTW
    delta_hat_NRTW  = 1-delta_hatnon_NRTW 
    theta_hat_NRTW = log(delta_hatnon_NRTW) + 2*(mu_hat_NRTW+sigmaSq_hat_NRTW)
    
    n_PTN = length(PTN)
    nonZero_PTN = PTN[which(PTN!=0)]
    nonzero_value_PTN = log(PTN[which(PTN!=0)])
    size_nonzero_PTN = length(nonZero_PTN)
    size_zero_PTN = n_PTN-length(nonZero_PTN)
    mu_hat_PTN = mean(nonzero_value_PTN)
    sigmaSq_hat_PTN  = var(nonzero_value_PTN)
    delta_hatnon_PTN = size_nonzero_PTN/n_PTN
    delta_hat_PTN  = 1-delta_hatnon_PTN 
    theta_hat_PTN = log(delta_hatnon_PTN) + 2*(mu_hat_PTN+sigmaSq_hat_PTN)
    
    
    n = data.frame()
    size_nonzero = data.frame()
    size_zero = data.frame()
    mu_hat = data.frame()
    sigmaSq_hat = data.frame()
    delta_hatnon = data.frame()
    delta_hat  = data.frame()
    theta_hat = data.frame()
    basicStat = data.frame()
    
    n = rbind(n_SK, n_SKSD, n_YL, n_NRTW, n_PTN)
    size_nonzero = rbind(size_nonzero_SK, size_nonzero_SKSD, size_nonzero_YL, size_nonzero_NRTW, size_nonzero_PTN)
    size_zero    = rbind(size_zero_SK, size_zero_SKSD, size_zero_YL, size_zero_NRTW, size_zero_PTN) 
    mu_hat     = rbind(mu_hat_SK, mu_hat_SKSD, mu_hat_YL, mu_hat_NRTW, mu_hat_PTN)
    sigmaSq_hat   = rbind(sigmaSq_hat_SK, sigmaSq_hat_SKSD, sigmaSq_hat_YL, sigmaSq_hat_NRTW, sigmaSq_hat_PTN)
    delta_hatnon = rbind(delta_hatnon_SK, delta_hatnon_SKSD, delta_hatnon_YL, delta_hatnon_NRTW, delta_hatnon_PTN)
    delta_hat    = 1-delta_hatnon
    theta_hat    = rbind(theta_hat_SK, theta_hat_SKSD, theta_hat_YL, theta_hat_NRTW, theta_hat_PTN)
    
    basicStat = cbind(n, size_nonzero, size_zero, mu_hat, sigmaSq_hat, delta_hat, theta_hat)
    colnames(basicStat) = c("n", "sizeNon", "sizeZero", "muHat", "sigmaSqHat", "deltaHat", "thetaHat")
    rownames(basicStat) = c("Songkhla", "Songkhla-Sadao", "Yala",  "Narathiwat", "Pattani")
    basicStat
    
    
    ##############################################################################
    # R.code to estimate AIC and BIC results                                     #
    ##############################################################################
    
    #SK
    SK.fit = SK[which(SK!=0)]
    norm.SK  = fitdist(SK.fit,"norm")
    summary.norm.SK = summary(norm.SK)
    AIC.norm.SK = summary.norm.SK$aic
    BIC.norm.SK = summary.norm.SK$bic
    
    lnorm.SK  = fitdist(SK.fit,"lnorm")
    summary.lnorm.SK = summary(lnorm.SK)
    AIC.lnorm.SK = summary.lnorm.SK$aic
    BIC.lnorm.SK = summary.lnorm.SK$bic
    
    logis.SK  = fitdist(SK.fit,"logis")
    summary.logis.SK = summary(logis.SK)
    AIC.logis.SK = summary.logis.SK$aic
    BIC.logis.SK = summary.logis.SK$bic
    
    expo.SK  = fitdist(SK.fit,"exp")
    summary.expo.SK = summary(expo.SK)
    AIC.expo.SK = summary.expo.SK$aic
    BIC.expo.SK = summary.expo.SK$bic
    
    cau.SK  = fitdist(SK.fit,"cauchy")
    summary.cau.SK = summary(cau.SK)
    AIC.cau.SK = summary.cau.SK$aic
    BIC.cau.SK = summary.cau.SK$bic
    
    
    #SKSD
    SKSD.fit = SKSD[which(SKSD!=0)]
    norm.SKSD  = fitdist(SKSD.fit,"norm")
    summary.norm.SKSD = summary(norm.SKSD)
    AIC.norm.SKSD = summary.norm.SKSD$aic
    BIC.norm.SKSD = summary.norm.SKSD$bic
    
    lnorm.SKSD  = fitdist(SKSD.fit,"lnorm")
    summary.lnorm.SKSD = summary(lnorm.SKSD)
    AIC.lnorm.SKSD = summary.lnorm.SKSD$aic
    BIC.lnorm.SKSD = summary.lnorm.SKSD$bic
    
    logis.SKSD  = fitdist(SKSD.fit,"logis")
    summary.logis.SKSD = summary(logis.SKSD)
    AIC.logis.SKSD = summary.logis.SKSD$aic
    BIC.logis.SKSD = summary.logis.SKSD$bic
    
    expo.SKSD  = fitdist(SKSD.fit,"exp")
    summary.expo.SKSD = summary(expo.SKSD)
    AIC.expo.SKSD = summary.expo.SKSD$aic
    BIC.expo.SKSD = summary.expo.SKSD$bic
    
    cau.SKSD  = fitdist(SKSD.fit,"cauchy")
    summary.cau.SKSD = summary(cau.SKSD)
    AIC.cau.SKSD = summary.cau.SKSD$aic
    BIC.cau.SKSD = summary.cau.SKSD$bic
    
    
    #YL
    YL.fit = YL[which(YL!=0)]
    norm.YL  = fitdist(YL.fit,"norm")
    summary.norm.YL = summary(norm.YL)
    AIC.norm.YL = summary.norm.YL$aic
    BIC.norm.YL = summary.norm.YL$bic
    
    lnorm.YL  = fitdist(YL.fit,"lnorm")
    summary.lnorm.YL = summary(lnorm.YL)
    AIC.lnorm.YL = summary.lnorm.YL$aic
    BIC.lnorm.YL = summary.lnorm.YL$bic
    
    logis.YL  = fitdist(YL.fit,"logis")
    summary.logis.YL = summary(logis.YL)
    AIC.logis.YL = summary.logis.YL$aic
    BIC.logis.YL = summary.logis.YL$bic
    
    expo.YL  = fitdist(YL.fit,"exp")
    summary.expo.YL = summary(expo.YL)
    AIC.expo.YL = summary.expo.YL$aic
    BIC.expo.YL = summary.expo.YL$bic
    
    cau.YL  = fitdist(YL.fit,"cauchy")
    summary.cau.YL = summary(cau.YL)
    AIC.cau.YL = summary.cau.YL$aic
    BIC.cau.YL = summary.cau.YL$bic
    
    
    
    #NRTW
    NRTW.fit = NRTW[which(NRTW!=0)]
    norm.NRTW  = fitdist(NRTW.fit,"norm")
    summary.norm.NRTW = summary(norm.NRTW)
    AIC.norm.NRTW = summary.norm.NRTW$aic
    BIC.norm.NRTW = summary.norm.NRTW$bic
    
    lnorm.NRTW  = fitdist(NRTW.fit,"lnorm")
    summary.lnorm.NRTW = summary(lnorm.NRTW)
    AIC.lnorm.NRTW = summary.lnorm.NRTW$aic
    BIC.lnorm.NRTW = summary.lnorm.NRTW$bic
    
    logis.NRTW  = fitdist(NRTW.fit,"logis")
    summary.logis.NRTW = summary(logis.NRTW)
    AIC.logis.NRTW = summary.logis.NRTW$aic
    BIC.logis.NRTW = summary.logis.NRTW$bic
    
    expo.NRTW  = fitdist(NRTW.fit,"exp")
    summary.expo.NRTW = summary(expo.NRTW)
    AIC.expo.NRTW = summary.expo.NRTW$aic
    BIC.expo.NRTW = summary.expo.NRTW$bic
    
    cau.NRTW  = fitdist(NRTW.fit,"cauchy")
    summary.cau.NRTW = summary(cau.NRTW)
    AIC.cau.NRTW = summary.cau.NRTW$aic
    BIC.cau.NRTW = summary.cau.NRTW$bic
    
    
    #PTN
    PTN.fit = PTN[which(PTN!=0)]
    norm.PTN  = fitdist(PTN.fit,"norm")
    summary.norm.PTN = summary(norm.PTN)
    AIC.norm.PTN = summary.norm.PTN$aic
    BIC.norm.PTN = summary.norm.PTN$bic
    
    lnorm.PTN  = fitdist(PTN.fit,"lnorm")
    summary.lnorm.PTN = summary(lnorm.PTN)
    AIC.lnorm.PTN = summary.lnorm.PTN$aic
    BIC.lnorm.PTN = summary.lnorm.PTN$bic
    
    logis.PTN  = fitdist(PTN.fit,"logis")
    summary.logis.PTN = summary(logis.PTN)
    AIC.logis.PTN = summary.logis.PTN$aic
    BIC.logis.PTN = summary.logis.PTN$bic
    
    expo.PTN  = fitdist(PTN.fit,"exp")
    summary.expo.PTN = summary(expo.PTN)
    AIC.expo.PTN = summary.expo.PTN$aic
    BIC.expo.PTN = summary.expo.PTN$bic
    
    cau.PTN  = fitdist(PTN.fit,"cauchy")
    summary.cau.PTN = summary(cau.PTN)
    AIC.cau.PTN = summary.cau.PTN$aic
    BIC.cau.PTN = summary.cau.PTN$bic
    
    
    AIC.SK = data.frame()
    AIC.SKSD = data.frame()
    AIC.YL = data.frame()
    AIC.NRTW = data.frame()
    AIC.PTN  = data.frame()
    AIC = data.frame()
    AIC.SK = rbind(AIC.norm.SK, AIC.lnorm.SK, AIC.logis.SK, AIC.expo.SK, AIC.cau.SK)
    AIC.SKSD = rbind(AIC.norm.SKSD, AIC.lnorm.SKSD, AIC.logis.SKSD, AIC.expo.SKSD, AIC.cau.SKSD)
    AIC.YL = rbind(AIC.norm.YL, AIC.lnorm.YL, AIC.logis.YL, AIC.expo.YL, AIC.cau.YL)
    AIC.NRTW = rbind(AIC.norm.NRTW, AIC.lnorm.NRTW, AIC.logis.NRTW, AIC.expo.NRTW, AIC.cau.NRTW)
    AIC.PTN = rbind(AIC.norm.PTN, AIC.lnorm.PTN, AIC.logis.PTN, AIC.expo.PTN, AIC.cau.PTN)
    AIC = cbind(AIC.SK, AIC.SKSD, AIC.YL, AIC.NRTW, AIC.PTN)
    
    rownames(AIC) = c("normal", "lognormal", "logistic",  "exponential", "cauchy")
    colnames(AIC) = c("Songkhla", "Songkhla-Sadao", "Yala",  "Narathiwat", "Pattani")
    AIC
    
    
    BIC.SK = data.frame()
    BIC.SKSD = data.frame()
    BIC.YL = data.frame()
    BIC.NRTW = data.frame()
    BIC.PTN  = data.frame()
    BIC = data.frame()
    BIC.SK = rbind(BIC.norm.SK, BIC.lnorm.SK, BIC.logis.SK, BIC.expo.SK, BIC.cau.SK)
    BIC.SKSD = rbind(BIC.norm.SKSD, BIC.lnorm.SKSD, BIC.logis.SKSD, BIC.expo.SKSD, BIC.cau.SKSD)
    BIC.YL = rbind(BIC.norm.YL, BIC.lnorm.YL, BIC.logis.YL, BIC.expo.YL, BIC.cau.YL)
    BIC.NRTW = rbind(BIC.norm.NRTW, BIC.lnorm.NRTW, BIC.logis.NRTW, BIC.expo.NRTW, BIC.cau.NRTW)
    BIC.PTN = rbind(BIC.norm.PTN, BIC.lnorm.PTN, BIC.logis.PTN, BIC.expo.PTN, BIC.cau.PTN)
    BIC = cbind(BIC.SK, BIC.SKSD, BIC.YL, BIC.NRTW, BIC.PTN)
    
    rownames(BIC) = c("normal", "lognormal", "logistic",  "exponential", "cauchy")
    colnames(BIC) = c("Songkhla", "Songkhla-Sadao", "Yala",  "Narathiwat", "Pattani")
    BIC