##################################################################################
# R.code for FGCI, LS, MOVER and highest posterior densities (HPDs)              #
##################################################################################


library(EnvStats)
library(base)
library(HDInterval)
library(LaplacesDemon)
library(invgamma)  
library(MASS)

# x    = random samples draw from delta-lognormal distribution (DLN)
# M    = number of simulation runs
# m    = number of generalized computation 
# n    = sample size 
# mu   = mean of normal distribution
# sigmaSq      = variance of normal distribution
# delta        = probability of additional zero
# alpha        = level of significance 
# Ex_LN        = mean of lognormal
# cv_LN        = coefficient of variation (CV) of lognormal
# Temp.size_nonzero = At least the sample size of nonzero = n-(n*delta)
# Var_theta        = variance of mean of DLN
# theta        = common mean of several DLNs   

# mu_FGPQ           = fiducial generalised pivotal quantity (FGPQ) of mu
# sigmaSq_FGPQ      = FGPQ of sigmaSq
# deltanon_FGPQ     = FGPQ of probability of non-zero values
# MeanDLN_FGPQ_i    = FGPQ of mean of DLN in the population i-th
# theta.FGCI        = FGPQ of common mean of DLNs
# L.FGCI            = lower limit of fiducial generalized confidence interval (FGCI)
# U.FGCI            = upper limit of FGCI
# CP.FGCI           = coverage probability of FGCI
# Length.FGCI       = average length of FGCI
# ACP.FGCI          = average coverage probability of FGCI      		                       
# Ex.length.FGCI    = average length of FGCI

# theta_hat   = the estimated common mean of DLNs
# var_thetahat.LS = the estimated variance of theta_hat
# L.LS        = lower limit of Method of Large sample (LS)
# U.LS        = upper limit of LS
# CP.LS       = coverage probability of LS
# Length.LS   = average length of LS
# ACP.LS      = average coverage probability of LS      		                       
# Ex.length.LS= average length of LS
# RAL.LS      = relative average length of LS

# L.mover           = lower limit of Method of variance estimates recovery (MOVER)
# U.mover           = upper limit of MOVER
# CP.mover          = coverage probability of MOVER
# Length.mover      = average length of MOVER
# ACP.mover         = average coverage probability of MOVER      		                       
# Ex.length.mover   = average length of MOVER
# RAL.mover         = relative average length of MOVER

# L.PB           = lower limit of parametric bootstrap (PB)
# U.PB           = upper limit of PB
# CP.PB          = coverage probability of PB
# Length.PB      = average length of PB
# ACP.PB         = average coverage probability of PB      		                       
# Ex.length.PB   = average length of PB
# RAL.PB         = relative average length of PB

# HPD.JR             = highest posterior density interval based on Jeffreys' Rule prior (HPD.JR)
# post.sigmaSq.JR    = posterior of sigmaSq based on Jeffreys' Rule prior
# post.mu.JR         = posterior of mu based on Jeffreys' Rule prior
# post.delta.JR      = posterior of delta based on Jeffreys' Rule prior
# postMeanDLN.JR     = posterior of delta-lognormal mean in the population i-th based on Jeffreys' Rule prior
# theta.JR           = posterior of common mean based on Jeffreys' Rule prior
# L.HPD.JR           = lower limit of HPD.JR
# U.HPD.JR           = upper limit of HPD.JR
# CP.HPD.JR          = coverage probability of HPD.JR
# Length.HPD.JR      = average length of HPD.JR
# ACP.HPD.JR         = average coverage probability of HPD.JR      		                       
# Ex.length.HPD.JR   = average length of HPD.JR

# HPD.NGB              = highest posterior density interval based on normal-gamma-beta prior (HPD-NGB)
# post.sigmaSq.NGB     = posterior of sigmaSq based on NGB prior
# post.mu.NGB          = posterior of mu based on NGB prior
# post.delta.NGB       = posterior of delta based on NGB prior 
# postMeanDLN.NGB      = posterior of delta-lognormal mean in the population i-th based on NGB prior
# theta.NGB            = posterior of common mean based on NGB prior
# L.HPD.NGB            = lower limit of HPD-NGB
# U.HPD.NGB            = upper limit of HPD-NGB
# CP.HPD.NGB           = coverage probability of HPD-NGB
# Length.HPD.NGB       = average length of HPD-NGB
# ACP.HPD.NGB          = average coverage probability of HPD-NGB      		                       
# Ex.length.HPD.NGB    = average length of HPD-NGB


##################################################################################
# k = 2 sample cases                                                             #
##################################################################################


#Function for generating random samples drawn from delta-lognormal distributions with k-th populations
xGen.k2 = function(k, nTotal, Ex_LN, cv_LN, deltaTotal){
  sum_nonzeroSq = rep(0,k)
  mu_hat = rep(0,k)
  sigmaSq_hat= rep(0,k)
  sigmaSq_hat.MLE= rep(0,k)
  delta_hatnon= rep(0,k)
  delta_hat= rep(0,k)
  size_non = rep(0,k)
  size_zero = rep(0,k)
  xlist = list()
  x_data = data.frame()
  
  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)])                         #find values is not zero
      Temp.size_nonzero =length(Temp.nonzero_value)                          #find size of nonzero
      if(Temp.size_nonzero < (nTotal[l]-nTotal[l]*deltaTotal[l])){
        next
      }
      nonzero_value = Temp.nonzero_value
      size_nonzero = Temp.size_nonzero
      sum_nonzeroSq[l] = sum(nonzero_value^2)
      mu_hat[l] = mean(nonzero_value)
      sigmaSq_hat[l] = var(nonzero_value)
      
      sigmaSq_hat.MLE[l] = (1/size_nonzero)*(sum(nonzero_value^2)-(size_nonzero*mu_hat[l]^2))
      
      delta_hatnon[l] = size_nonzero/nTotal[l]
      delta_hat[l] = 1- delta_hatnon[l]
      size_non[l] = size_nonzero
      size_zero[l] = nTotal[l]-size_nonzero
      
      tmp_x = data.frame(x)
      tmp_x$xlabel = l
      names(tmp_x) = c("data","xlabel")
      x_data = rbind(x_data,tmp_x)
      
      break
    }
  }
  codeN = paste(nTotal, collapse = '')
  filename = paste("result_x_K",k,"_",codeN,".csv", sep='')
  write.csv(x_data,filename,row.names = FALSE)
  return(cbind(size_non,size_zero, mu_hat, sigmaSq_hat, sigmaSq_hat.MLE, delta_hat, sum_nonzeroSq))	
}


#Function for computing CIs for the common mean of several delta-lognormal distributions
CI.commonMean.k2 = function(M=5000, m=2500, n=100, k=2, 
                         n1=30, n2=30, 
                         delta1=0.1, delta2=0.2, 
                         sigmaSq1=1, sigmaSq2=2, alpha=0.05){
  
  starttime = proc.time()
  
  nTotal = c(n1,n2)
  sigmaSqTotal = c(sigmaSq1,sigmaSq2)
  meanTotal =  c(0,0)
  deltaTotal = c(delta1,delta2)
  size_nonzeroTotal = nTotal*(1-deltaTotal)
  k = length(nTotal)
  
  Ex_LN = exp(meanTotal+(sigmaSqTotal/2))
  cv_LN = sqrt(exp(sigmaSqTotal)-1)
  
  #Variance of mean of DLN
  var1 = exp((2*meanTotal)+sigmaSqTotal)
  var2 = deltaTotal*(1-deltaTotal)
  var3 = 0.5*(1-deltaTotal)*((2*sigmaSqTotal)+sigmaSqTotal^2)
  Var_theta = var1*(var2+var3)/nTotal
  
  #The common mean of several DLNs
  w_i = 1/Var_theta
  theta = sum(w_i*(1-deltaTotal)*exp(meanTotal+(sigmaSqTotal/2)))/sum(w_i)
  

  #FGCI
  theta.FGCI = matrix(rep(0,m))
  L.FGCI = numeric()
  U.FGCI = numeric()
  CP.FGCI = rep(0,M)
  Length.FGCI = rep(0,M)
  
  #LS
  L.LS = numeric()
  U.LS  = numeric()
  CP.LS  = rep(0,M) 
  Length.LS  = rep(0,M)
  
  #MOVER
  L.mover= numeric()
  U.mover= numeric()
  CP.mover= rep(0,M)
  Length.mover= rep(0,M)
  
  #PB
  mu_hat.PB.i = rep(0,k)
  sigmaSq_hat.PB.i = rep(0,k)
  delta_hat.PB.i = rep(0,k)
  size_nonzero.PB.i = rep(0,k)
  Q.PB = rep(0,m)
  L.PB = numeric()
  U.PB = numeric()
  CP.PB = rep(0,M)
  Length.PB = rep(0,M)
  
  #HPD.JR
  post.mu.JR = rep(0,k)
  post.sigmaSq.JR = rep(0,k)
  post.delta.JR = rep(0,k)
  postMeanDLN.JR = rep(0,k)
  theta.JR = rep(0,m)
  L.HPD.JR = numeric()
  U.HPD.JR = numeric()
  CP.HPD.JR = rep(0,M)
  Length.HPD.JR = rep(0,M)
  
  #HPD.NGB
  post.mu.NGB = rep(0,k)
  post.sigmaSq.NGB = rep(0,k)
  post.delta.NGB = rep(0,k)
  postMeanDLN.NGB = rep(0,k) 
  theta.NGB = rep(0,m)
  L.HPD.NGB = numeric()
  U.HPD.NGB = numeric()
  CP.HPD.NGB = rep(0,M)
  Length.HPD.NGB = rep(0,M)
  
  i = 0
  while(i < M){                                                   # begin loop i
    
    
    data = xGen.k2(k, nTotal, Ex_LN, cv_LN,deltaTotal)
    codeN = paste(nTotal, collapse = '')
    filename = paste("result_x_K",k,"_",codeN,".csv", sep='')
    x_data = read.csv(filename)
    
    
    size_nonzero = data[,1]
    size_zero = data[,2]
    mu_hat = data[,3]
    sigmaSq_hat  = data[,4]
    sigmaSq_hat.MLE  = data[,5]
    delta_hat= data[,6]
    sum_nonzeroSq= data[,7]
    n = nTotal
    
    if(sum(delta_hat==0)>0){next}
    
    i=i+1
    if(i%%1000==0){
      print(i)
    }
    
    ##############################################################################
    #                         Algorithm 1: FGCI                                   #
    ##############################################################################
    e = 1
    while(e <= m){
      
        USq = rchisq(k,size_nonzero-1)/(size_nonzero-1)
        mu_FGPQ = mu_hat - ((rnorm(k)/sqrt(USq))*sqrt(sigmaSq_hat/size_nonzero))
        sigmaSq_FGPQ = sigmaSq_hat/USq
        deltanon_FGPQ = rbeta(k,size_nonzero+0.5, size_zero+0.5)
        delta_FGPQ = 1-deltanon_FGPQ
        MeanDLN_FGPQ_i  = deltanon_FGPQ*exp(mu_FGPQ+(sigmaSq_FGPQ/2))
      
        #Variance of meanDLN_hat
        var1.FGCI = exp((2*mu_FGPQ)+sigmaSq_FGPQ)
        var2.FGCI = delta_FGPQ*(1-delta_FGPQ)
        var3.FGCI = 0.5*(1-delta_FGPQ)*((2*sigmaSq_FGPQ)+sigmaSq_FGPQ^2)
        Var_thetahat.FGCI = var1.FGCI*(var2.FGCI+var3.FGCI)/n
      
        w_i.FGCI = 1/Var_thetahat.FGCI
        theta.FGCI[e] = sum(w_i.FGCI*MeanDLN_FGPQ_i)/sum(w_i.FGCI)
        e=e+1
    }
    
    L.FGCI = quantile(theta.FGCI,0.025)                                       
    U.FGCI = quantile(theta.FGCI,0.975)      
    CP.FGCI[i] = ifelse(L.FGCI< theta && theta <U.FGCI,1,0)
    Length.FGCI[i] = U.FGCI - L.FGCI
    
    #####################################################################
    #                     Algorithm 2: LS                               #
    #####################################################################
    
    k_df =size_nonzero-1
    t = sigmaSq_hat/2
    Bessel1 = 1+((size_nonzero-1)/size_nonzero)*t
    Bessel2 = ((size_nonzero-1)^3/(factorial(2)*size_nonzero^2))*(t^2/(size_nonzero+1))
    Bessel3 = ((size_nonzero-1)^5/(factorial(3)*size_nonzero^3))*(t^3/((size_nonzero+1)*(size_nonzero+3)))
    Bessel4 = ((size_nonzero-1)^7/(factorial(4)*size_nonzero^4))*(t^4/((size_nonzero+1)*(size_nonzero+3)*(size_nonzero+5)))
    Bessel5 = ((size_nonzero-1)^9/(factorial(5)*size_nonzero^5))*(t^5/((size_nonzero+1)*(size_nonzero+3)*(size_nonzero+5)*(size_nonzero+7)))
    Bessel= Bessel1+Bessel2+Bessel3+Bessel4+Bessel5
    theta_hat_i = (1-delta_hat)*exp(mu_hat)*Bessel
    
    #Variance of meanDLN_hat
    var1.LS  = exp((2*mu_hat)+sigmaSq_hat)
    var2.LS  = delta_hat*(1-delta_hat)
    var3.LS = 0.5*(1-delta_hat)*((2*sigmaSq_hat)+sigmaSq_hat^2)
    var_thetahat.LS = var1.LS *(var2.LS +var3.LS)/n
    
    w_i = 1/var_thetahat.LS
    theta_hat = sum(w_i*theta_hat_i)/sum(w_i)
    z =  qnorm(1-(alpha/2),0,1)
    sd.LS = sqrt(1/(sum(w_i)))
    
    L.LS = theta_hat-(z*sd.LS)                                      
    U.LS = theta_hat+(z*sd.LS)           
    CP.LS[i] = ifelse(L.LS< theta && theta<U.LS,1,0)                 
    Length.LS[i] = U.LS-L.LS                                        
    
    ##############################################################################
    #                          Algorithm 3:  MOVER                               #
    ##############################################################################
    
    #Wilson interval for delta
    r = qnorm(1-(alpha/2))
    deltanon_hat = 1 - delta_hat
    CI.delta1 = deltanon_hat +(0.5*r^2/n)
    CI.delta2 = sqrt((deltanon_hat*(1-deltanon_hat)+ (0.25*r^2/n))/n)
    CI.delta3 = 1+(r^2/n) 
    l.delta = log((CI.delta1-(r*CI.delta2))/CI.delta3)
    u.delta = log((CI.delta1+(r*CI.delta2))/CI.delta3)
   

    #CI for sigmaSq
    df = size_nonzero-1
    locate_L = qchisq(alpha/2,df) 
    locate_U = qchisq(1-(alpha/2),df) 
    
    #MOVER CI for mu + sigmaSq/2
    sd_mu = qnorm(1-(alpha/2))*sqrt(sigmaSq_hat/size_nonzero)
    sd_sigmaSq_L = (sigmaSq_hat/2)*(1-((size_nonzero-1)/(locate_U)))
    sd_sigmaSq_U = (sigmaSq_hat/2)*(((size_nonzero-1)/(locate_L))-1)
    L_the2_3= (mu_hat+sigmaSq_hat/2)-sqrt(sd_mu^2 + sd_sigmaSq_L^2)
    U_the2_3= (mu_hat+sigmaSq_hat/2)+sqrt(sd_mu^2 + sd_sigmaSq_U^2)
    
    #MOVER CI for log(1-delta_hat) + mu_hat + sigmaSq_hat/2
    l.MOVER = exp((log(1-delta_hat)+(mu_hat + sigmaSq_hat/2))-sqrt((log(1-delta_hat)-l.delta)^2+((mu_hat + sigmaSq_hat/2)-L_the2_3)^2))
    u.MOVER = exp((log(1-delta_hat)+(mu_hat + sigmaSq_hat/2))+sqrt((u.delta-log(1-delta_hat))^2+(U_the2_3-(mu_hat + sigmaSq_hat/2))^2))

   #Variance of estimated mean of DLN
    var1.MOVER  = exp((2*mu_hat)+sigmaSq_hat)
    var2.MOVER  = delta_hat*(1-delta_hat)
    var3.MOVER = 0.5*(1-delta_hat)*((2*sigmaSq_hat)+sigmaSq_hat^2)
    var_thetahat.MOVER = var1.MOVER *(var2.MOVER +var3.MOVER)/n
    w_i = 1/var_thetahat.MOVER
    
    #The estimated mean of DLN
    t = sigmaSq_hat/2
    Bessel1 = 1+((size_nonzero-1)/size_nonzero)*t
    Bessel2 = ((size_nonzero-1)^3/(factorial(2)*size_nonzero^2))*(t^2/(size_nonzero+1))
    Bessel3 = ((size_nonzero-1)^5/(factorial(3)*size_nonzero^3))*(t^3/((size_nonzero+1)*(size_nonzero+3)))
    Bessel4 = ((size_nonzero-1)^7/(factorial(4)*size_nonzero^4))*(t^4/((size_nonzero+1)*(size_nonzero+3)*(size_nonzero+5)))
    Bessel5 = ((size_nonzero-1)^9/(factorial(5)*size_nonzero^5))*(t^5/((size_nonzero+1)*(size_nonzero+3)*(size_nonzero+5)*(size_nonzero+7)))
    Bessel= Bessel1+Bessel2+Bessel3+Bessel4+Bessel5
    theta_i.MOVER = (1-delta_hat)*exp(mu_hat)*Bessel
    
    #MOVER CI for the common mean of several DLNs
    Temp.l.MOVER = (sum(w_i*theta_i.MOVER)/sum(w_i))- sqrt(sum(w_i^2*(theta_i.MOVER-l.MOVER)^2)/sum(w_i^2))       #compute lower based on MOVER VST
    Temp.u.MOVER = (sum(w_i*theta_i.MOVER)/sum(w_i))+ sqrt(sum(w_i^2*(theta_i.MOVER-u.MOVER)^2)/sum(w_i^2))      #compute upper based on MOVER VST
    
    l.mover= ifelse(w_i>0, Temp.l.MOVER, Temp.u.MOVER)
    u.mover= ifelse(w_i>0, Temp.u.MOVER, Temp.l.MOVER)
    
    L.mover= (sum(w_i*theta_i.MOVER)/sum(w_i))- sqrt(sum(w_i^2*(theta_i.MOVER-l.mover)^2)/sum(w_i^2))       #compute lower based on MOVER-Wilson
    U.mover= (sum(w_i*theta_i.MOVER)/sum(w_i))+ sqrt(sum(w_i^2*(theta_i.MOVER-u.mover)^2)/sum(w_i^2))      #compute upper based on MOVER-Wilson
    CP.mover[i] = ifelse(L.mover< theta && theta< U.mover,1,0)                  #compute probability
    Length.mover[i] = U.mover- L.mover                                         #compute length based on MOVER-Wilson
    
    #####################################################################
    #                  Algorithm 4: PB                                  #
    #####################################################################
    
    #The estimate SigmaSq
    w_hat.i = size_nonzero/sigmaSq_hat.MLE
    lnMeanDLN_hat = sum(w_hat.i*(mu_hat+(sigmaSq_hat.MLE/2)+log(1-delta_hat)))/sum(w_hat.i)
    AA = mu_hat - (lnMeanDLN_hat-log(1-delta_hat))
    
    #The MLE of SigmaSq
    sigmaSq_hat.MLEPB = 2*(sqrt(1 + sigmaSq_hat.MLE + AA^2)-1)
    w_hat.i.MLE = size_nonzero/sigmaSq_hat.MLEPB
    delta_hat.i.MLE = size_zero/(n+w_hat.i.MLE)
    lnMeanDLN_hat.MLE = sum(w_hat.i.MLE*(mu_hat+(sigmaSq_hat.MLEPB/2)+log(1-delta_hat.i.MLE)))/sum(w_hat.i.MLE)  
    
    h = 1
    while(h <= m){
      for(i_k in 1:k){
        x_DLN = x_data[which(x_data$xlabel==i_k),]
        
        A = sample(1:n[i_k],size=n[i_k],replace=TRUE)
        x_DLN.BS = x_DLN[A,]
        
        nonzero_value.BS = log(x_DLN.BS[which(x_DLN.BS$data!=0),'data'])
        size_nonzero.BS =length(nonzero_value.BS)
        size_zero.BS = n[i_k]-size_nonzero.BS
        mu_hat.BS = mean(nonzero_value.BS)
        sigmaSq_hat.BS = var(nonzero_value.BS)
        deltahat.BS = size_zero.BS/n[i_k]
        
        mu_hat.PB.i[i_k] = rnorm(1, mu_hat.BS, sqrt(sigmaSq_hat.BS/size_nonzero.BS))
        sigmaSq_hat.PB.i[i_k] = sigmaSq_hat.BS*rchisq(1,size_nonzero.BS-1)/(size_nonzero.BS-1)
        delta_hat.PB.i[i_k] = rbeta(1, size_zero.BS+0.5, size_nonzero.BS+0.5)
        size_nonzero.PB.i[i_k] = n[i_k]*(1-delta_hat.PB.i[i_k])
  
      }
      
      w.PB.i = size_nonzero.PB.i/sigmaSq_hat.PB.i
      lnMeanDLN.PB = sum(w.PB.i*(mu_hat.PB.i+(sigmaSq_hat.PB.i/2)+log(1-delta_hat.PB.i)))/sum(w.PB.i)
      
      AA.MLE = mu_hat.PB.i - (lnMeanDLN.PB-log(1-delta_hat.PB.i))
      sigmaSq_hat.MLE.PB = 2*(sqrt(1 + sigmaSq_hat.PB.i + AA.MLE^2)-1)
      w_hat.i.MLE.PB = size_nonzero.PB.i/sigmaSq_hat.MLE.PB
      delta_hat.i.MLE.PB = (n-size_nonzero.PB.i)/(n+w_hat.i.MLE.PB)
      lnMeanDLN_hat.MLE.PB = sum(w_hat.i.MLE.PB*(mu_hat.PB.i+(sigmaSq_hat.MLE.PB/2)+log(1-delta_hat.i.MLE.PB)))/sum(w_hat.i.MLE.PB)  
      
      Q.PB[h] = (lnMeanDLN_hat.MLE.PB-lnMeanDLN_hat.MLE)^2*sum(w_hat.i.MLE.PB)
      if(is.na(Q.PB[h])){next}
      h = h+1
    }
    
    q.PB = quantile(Q.PB, 1-alpha)
    L.PB = exp(lnMeanDLN_hat.MLE - sqrt(q.PB/sum(w_hat.i.MLE)))
    U.PB = exp(lnMeanDLN_hat.MLE + sqrt(q.PB/sum(w_hat.i.MLE)))
    CP.PB[i] = ifelse(L.PB < theta && theta < U.PB,1,0)
    Length.PB[i]= U.PB - L.PB
    
    
    #####################################################################
    #               Algorithm 5: HPD-JR and HPD-NGB                     #
    #####################################################################
    
    ############################### HPD-JR ##############################
    r = 0
    while(r < m){
      
      v = size_nonzero+1
      shap = v/2
      scal = v*sigmaSq_hat/2
      post.sigmaSq.JR = rinvgamma(k, shap, scal)
      
      Tem.post.mu.JR = numeric()
      for(k in 1:length(post.sigmaSq.JR)){
        Tem.post.mu.JR[k] = rnorm(1, mu_hat, sqrt(post.sigmaSq.JR[k]/size_nonzero))
      }
      post.mu.JR = Tem.post.mu.JR
      post.delta.JR  = rbeta(k, size_zero+0.5, size_nonzero+1.5)   
      postMeanDLN.JR = (1-post.delta.JR)*exp(post.mu.JR + post.sigmaSq.JR/2)
      
      var1.JR = exp((2*post.mu.JR)+post.sigmaSq.JR)
      var2.JR = post.delta.JR*(1-post.delta.JR)
      var3.JR = 0.5*(1-post.delta.JR)*((2*post.sigmaSq.JR)+post.sigmaSq.JR^2)
      var_thetahat.JR = var1.JR*(var2.JR+var3.JR)/n
      
      w_i.JR = 1/var_thetahat.JR
      theta.JR [r] = sum(w_i.JR*postMeanDLN.JR)/sum(w_i.JR)
      
      r=r+1
    }
    
    HPD.JR = hdi(theta.JR, credMass = 1-alpha)
    L.HPD.JR  = quantile(HPD.JR ,0.025)                                       
    U.HPD.JR  = quantile(HPD.JR ,0.975)           
    CP.HPD.JR [i] = ifelse(L.HPD.JR < theta && theta < U.HPD.JR ,1,0)
    Length.HPD.JR [i] = U.HPD.JR -L.HPD.JR     
    
    
    ############################### HPD-NGB ##############################
    
    s = 0
    while(s < m){
      
      alpha_n = (size_nonzero-1)/2
      beta_n = (sum_nonzeroSq-(size_nonzero*mu_hat^2))/2
      post.sigmaSq.NGB = rinvgamma(k, alpha_n, beta_n)
      post.mu.NGB = rst(k, mu_hat, sqrt(sigmaSq_hat/size_nonzero), v)
      post.delta.NGB  = rbeta(k, size_zero+0.5, size_nonzero+0.5)   
      postMeanDLN.NGB  = (1-post.delta.NGB)*exp(post.mu.NGB + post.sigmaSq.NGB/2)
      
      var1.NGB = exp((2*post.mu.NGB)+post.sigmaSq.NGB)
      var2.NGB = post.delta.NGB*(1-post.delta.NGB)
      var3.NGB = 0.5*(1-post.delta.NGB)*((2*post.sigmaSq.NGB)+post.sigmaSq.NGB^2)
      var_thetahat.NGB = var1.NGB*(var2.NGB+var3.NGB)/n
      
      w_i.NGB = 1/var_thetahat.NGB
      theta.NGB[s] = sum(w_i.NGB*postMeanDLN.NGB)/sum(w_i.NGB)
      s=s+1
    }
    
    HPD.NGB = hdi(theta.NGB, credMass = 1-alpha)
    L.HPD.NGB  = quantile(HPD.NGB ,0.025)                                       
    U.HPD.NGB  = quantile(HPD.NGB ,0.975)           
    CP.HPD.NGB[i] = ifelse(L.HPD.NGB < theta && theta < U.HPD.NGB ,1,0)
    Length.HPD.NGB[i] = U.HPD.NGB -L.HPD.NGB     
    
    
  }                                                                                #end loop i
 
  ACP.FGCI = mean(CP.FGCI)	                                                    
  Ex.length.FGCI = mean(Length.FGCI)                                          
  
  ACP.LS = mean(CP.LS)	                                                          
  Ex.length.LS = mean(Length.LS)                                                
  
  ACP.mover= mean(CP.mover)	                                               
  Ex.length.mover= mean(Length.mover)                                       
  
  ACP.PB = mean(CP.PB)	                                               
  Ex.length.PB = mean(Length.PB)                                      
  
  ACP.HPD.JR = mean(CP.HPD.JR)	                                               
  Ex.length.HPD.JR = mean(Length.HPD.JR)                                       
  
  ACP.HPD.NGB = mean(CP.HPD.NGB)	                                               
  Ex.length.HPD.NGB = mean(Length.HPD.NGB)                                      
  
  
  elapseTime = proc.time()-starttime                                   
  
  cat("k =",k,"M =",M,"m =",m, "nTotal =",nTotal, "meanTotal=",meanTotal, "sigmaSqTotal=",sigmaSqTotal,"deltaTotal =",deltaTotal, "theta=",theta,"\n",
     
      "CP|FGCI. =",ACP.FGCI ,"Length|FGCI  =",Ex.length.FGCI,"\n",
      "CP|LS =",ACP.LS,"Length|LS =",Ex.length.LS,"\n",
      "CP|mover =",ACP.mover,"Length|mover =",Ex.length.mover,"\n",
      "CP|PB =",ACP.PB,"Length|PB =",Ex.length.PB,"\n",
      "CP|HPD.JR =",ACP.HPD.JR,"Length|HPD.JR =",Ex.length.HPD.JR,"\n",
      "CP|HPD.NGB =",ACP.HPD.NGB,"Length|HPD.NGB =",Ex.length.HPD.NGB)
  
  
  result = c("k =",k,"M =",M,"m =",m, "nTotal =",nTotal, "meanTotal=",meanTotal, "sigmaSqTotal=",sigmaSqTotal,"deltaTotal =",deltaTotal, "theta=",theta,
             
            
             "CP|FGCI =",ACP.FGCI,
             "CP|LS =",ACP.LS,
             "CP|mover =",ACP.mover,
             "CP|PB =",ACP.PB,
             "CP|HPD.JR =",ACP.HPD.JR,
             "CP|HPD.NGB =",ACP.HPD.NGB,
             
             "Length|FGCI  =",Ex.length.FGCI,
             "Length|LS =",Ex.length.LS,
             "Length|mover =",Ex.length.mover,
             "Length|PB =",Ex.length.PB,
             "Length|HPD.JR =",Ex.length.HPD.JR,
             "Length|HPD.NGB =",Ex.length.HPD.NGB,
             
             "Time =",elapseTime["elapsed"])
  return(result)
  
  
}



##################################################################################
# k = 5 sample cases                                                             #
##################################################################################

xGen.k5 = function(k,nTotal,Ex_LN,cv_LN,deltaTotal){
  sum_nonzeroSq = rep(0,k)
  mu_hat = rep(0,k)
  sigmaSq_hat= rep(0,k)
  sigmaSq_hat.MLE= rep(0,k)
  delta_hatnon= rep(0,k)
  delta_hat= rep(0,k)
  size_non = rep(0,k)
  size_zero = rep(0,k)
  xlist = list()
  x_data = data.frame()
  
  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)])                         #find values is not zero
      Temp.size_nonzero =length(Temp.nonzero_value)                          #find size of nonzero
      if(Temp.size_nonzero < (nTotal[l]-nTotal[l]*deltaTotal[l])){
        next
      }
      nonzero_value = Temp.nonzero_value
      size_nonzero = Temp.size_nonzero
      sum_nonzeroSq[l] = sum(nonzero_value^2)
      mu_hat[l] = mean(nonzero_value)
      sigmaSq_hat[l] = var(nonzero_value)
      
      sigmaSq_hat.MLE[l] = (1/size_nonzero)*(sum(nonzero_value^2)-(size_nonzero*mu_hat[l]^2))
      
      delta_hatnon[l] = size_nonzero/nTotal[l]
      delta_hat[l] = 1- delta_hatnon[l]
      size_non[l] = size_nonzero
      size_zero[l] = nTotal[l]-size_nonzero
      
      tmp_x = data.frame(x)
      tmp_x$xlabel = l
      names(tmp_x) = c("data","xlabel")
      x_data = rbind(x_data,tmp_x)
      
      break
    }
  }
  codeN = paste(nTotal, collapse = '')
  filename = paste("result_x_K",k,"_",codeN,".csv", sep='')
  write.csv(x_data,filename,row.names = FALSE)
  return(cbind(size_non,size_zero, mu_hat, sigmaSq_hat, sigmaSq_hat.MLE, delta_hat, sum_nonzeroSq))	
}

CI.CommonMean.k5 = function(M=5000,m=2500,k=5, 
                            n1=30, n2=30, n3=30, n4=30, n5=30, 
                            delta1=0.05, delta2=0.1, delta3=0.1, delta4=0.2, delta5=0.2, 
                            sigmaSq1=1, sigmaSq2=1, sigmaSq3=2, sigmaSq4=2, sigmaSq5=2,
                            alpha=0.05){
  starttime = proc.time()
  
  
  nTotal = c(n1,n2,n3,n4,n5)
  sigmaSqTotal = c(sigmaSq1,sigmaSq2,sigmaSq3,sigmaSq4,sigmaSq5)
  meanTotal =  c(0,0,0,0,0)
  deltaTotal = c(delta1,delta2,delta3,delta4,delta5)
  size_nonzeroTotal = nTotal*(1-deltaTotal)
  k = length(nTotal)
  
  Ex_LN = exp(meanTotal+(sigmaSqTotal/2))
  cv_LN = sqrt(exp(sigmaSqTotal)-1)
  
  var1 = exp((2*meanTotal)+sigmaSqTotal)
  var2 = deltaTotal*(1-deltaTotal)
  var3 = 0.5*(1-deltaTotal)*((2*sigmaSqTotal)+sigmaSqTotal^2)
  Var_theta = var1*(var2+var3)/nTotal
  
  w_i = 1/Var_theta
  theta = sum(w_i*(1-deltaTotal)*exp(meanTotal+(sigmaSqTotal/2)))/sum(w_i)
  
  #FGCI
  theta.FGCI = matrix(rep(0,m))
  L.FGCI = numeric()
  U.FGCI = numeric()
  CP.FGCI = rep(0,M)
  Length.FGCI = rep(0,M)
  
  #LS
  L.LS = numeric()
  U.LS  = numeric()
  CP.LS  = rep(0,M) 
  Length.LS  = rep(0,M)
  
  #MOVER
  L.mover = numeric()
  U.mover = numeric()
  CP.mover = rep(0,M)
  Length.mover = rep(0,M)
  
  #PB
  mu_hat.PB.i = rep(0,k)
  sigmaSq_hat.PB.i = rep(0,k)
  delta_hat.PB.i = rep(0,k)
  size_nonzero.PB.i = rep(0,k)
  Q.PB = rep(0,m)
  L.PB = numeric()
  U.PB = numeric()
  CP.PB = rep(0,M)
  Length.PB = rep(0,M)
  
  #HPD-JR
  post.mu.JR = rep(0,k)
  post.sigmaSq.JR = rep(0,k)
  post.delta.JR = rep(0,k)
  postMeanDLN.JR = rep(0,k)
  theta.JR = rep(0,m)
  L.HPD.JR = numeric()
  U.HPD.JR = numeric()
  CP.HPD.JR = rep(0,M)
  Length.HPD.JR = rep(0,M)
  
  #HPD-NGB
  post.mu.NGB = rep(0,k)
  post.sigmaSq.NGB = rep(0,k)
  post.delta.NGB = rep(0,k)
  postMeanDLN.NGB = rep(0,k) 
  theta.NGB = rep(0,m)
  L.HPD.NGB = numeric()
  U.HPD.NGB = numeric()
  CP.HPD.NGB = rep(0,M)
  Length.HPD.NGB = rep(0,M)
  
  
  i = 0
  while(i < M){                                                   # begin loop i
    
    
    data = xGen.k5(k, nTotal, Ex_LN, cv_LN,deltaTotal)
    codeN = paste(nTotal, collapse = '')
    filename = paste("result_x_K",k,"_",codeN,".csv", sep='')
    x_data = read.csv(filename)
    
    
    size_nonzero = data[,1]
    size_zero = data[,2]
    mu_hat = data[,3]
    sigmaSq_hat  = data[,4]
    sigmaSq_hat.MLE  = data[,5]
    delta_hat= data[,6]
    sum_nonzeroSq= data[,7]
    n = nTotal
    
    if(sum(delta_hat==0)>0){next}
    i=i+1
    if(i%%1000==0){
      print(i)
    }
    
    
    ##############################################################################
    #                         Algorithm 1: FGCI                                  #
    ##############################################################################
    e = 1
    while(e <= m){
      
      USq = rchisq(k,size_nonzero-1)/(size_nonzero-1)
      mu_FGPQ = mu_hat - ((rnorm(k)/sqrt(USq))*sqrt(sigmaSq_hat/size_nonzero))
      sigmaSq_FGPQ = sigmaSq_hat/USq
      deltanon_FGPQ = rbeta(k,size_nonzero+0.5, size_zero+0.5)
      delta_FGPQ = 1-deltanon_FGPQ
      MeanDLN_FGPQ_i  = deltanon_FGPQ*exp(mu_FGPQ+(sigmaSq_FGPQ/2))
      
      #Variance of meanDLN_hat based on Owen and DeRouen 1980
      var1.FGCI = exp((2*mu_FGPQ)+sigmaSq_FGPQ)
      var2.FGCI = delta_FGPQ*(1-delta_FGPQ)
      var3.FGCI = 0.5*(1-delta_FGPQ)*((2*sigmaSq_FGPQ)+sigmaSq_FGPQ^2)
      Var_thetahat.FGCI = var1.FGCI*(var2.FGCI+var3.FGCI)/n
      
      w_i.FGCI = 1/Var_thetahat.FGCI
      theta.FGCI[e] = sum(w_i.FGCI*MeanDLN_FGPQ_i)/sum(w_i.FGCI)
      e=e+1
    }
    
    L.FGCI = quantile(theta.FGCI,0.025)                                       
    U.FGCI = quantile(theta.FGCI,0.975)      
    CP.FGCI[i] = ifelse(L.FGCI< theta && theta <U.FGCI,1,0)
    Length.FGCI[i] = U.FGCI - L.FGCI
    
    #####################################################################
    #                     Algorithm 2: LS                               #
    #####################################################################
    
    k_df =size_nonzero-1
    
    #Bessel function
    t = sigmaSq_hat/2
    
    Bessel1 = 1+((size_nonzero-1)/size_nonzero)*t
    Bessel2 = ((size_nonzero-1)^3/(factorial(2)*size_nonzero^2))*(t^2/(size_nonzero+1))
    Bessel3 = ((size_nonzero-1)^5/(factorial(3)*size_nonzero^3))*(t^3/((size_nonzero+1)*(size_nonzero+3)))
    Bessel4 = ((size_nonzero-1)^7/(factorial(4)*size_nonzero^4))*(t^4/((size_nonzero+1)*(size_nonzero+3)*(size_nonzero+5)))
    Bessel5 = ((size_nonzero-1)^9/(factorial(5)*size_nonzero^5))*(t^5/((size_nonzero+1)*(size_nonzero+3)*(size_nonzero+5)*(size_nonzero+7)))
    Bessel= Bessel1+Bessel2+Bessel3+Bessel4+Bessel5
    theta_hat_i = (1-delta_hat)*exp(mu_hat)*Bessel
    
    var1.LS  = exp((2*mu_hat)+sigmaSq_hat)
    var2.LS  = delta_hat*(1-delta_hat)
    var3.LS = 0.5*(1-delta_hat)*((2*sigmaSq_hat)+sigmaSq_hat^2)
    var_thetahat.LS = var1.LS *(var2.LS +var3.LS)/n
    
    w_i = 1/var_thetahat.LS
    theta_hat = sum(w_i*theta_hat_i)/sum(w_i)
    z =  qnorm(1-(alpha/2),0,1)
    sd.LS = sqrt(1/(sum(w_i)))
    
    L.LS = theta_hat-(z*sd.LS)                                      
    U.LS = theta_hat+(z*sd.LS)           
    CP.LS[i] = ifelse(L.LS< theta && theta<U.LS,1,0)                 
    Length.LS[i] = U.LS-L.LS                                        
    
    ##############################################################################
    #                          Algorithm 3:  MOVER                               #
    ##############################################################################
    
    #Wilson interval for delta
    r = qnorm(1-(alpha/2))
    deltanon_hat = 1 - delta_hat
    CI.delta1 = deltanon_hat +(0.5*r^2/n)
    CI.delta2 = sqrt((deltanon_hat*(1-deltanon_hat)+ (0.25*r^2/n))/n)
    CI.delta3 = 1+(r^2/n) 
    l.delta = log((CI.delta1-(r*CI.delta2))/CI.delta3)
    u.delta = log((CI.delta1+(r*CI.delta2))/CI.delta3)
    
    
    #CI for sigmaSq
    df = size_nonzero-1
    locate_L = qchisq(alpha/2,df) 
    locate_U = qchisq(1-(alpha/2),df) 
    
    #MOVER CI for mu + sigmaSq/2
    sd_mu = qnorm(1-(alpha/2))*sqrt(sigmaSq_hat/size_nonzero)
    sd_sigmaSq_L = (sigmaSq_hat/2)*(1-((size_nonzero-1)/(locate_U)))
    sd_sigmaSq_U = (sigmaSq_hat/2)*(((size_nonzero-1)/(locate_L))-1)
    L_the2_3= (mu_hat+sigmaSq_hat/2)-sqrt(sd_mu^2 + sd_sigmaSq_L^2)
    U_the2_3= (mu_hat+sigmaSq_hat/2)+sqrt(sd_mu^2 + sd_sigmaSq_U^2)
    
    #MOVER CI for log(1-delta_hat) + mu_hat + sigmaSq_hat/2
    l.MOVER = exp((log(1-delta_hat)+(mu_hat + sigmaSq_hat/2))-sqrt((log(1-delta_hat)-l.delta)^2+((mu_hat + sigmaSq_hat/2)-L_the2_3)^2))
    u.MOVER = exp((log(1-delta_hat)+(mu_hat + sigmaSq_hat/2))+sqrt((u.delta-log(1-delta_hat))^2+(U_the2_3-(mu_hat + sigmaSq_hat/2))^2))
    
    #Variance of estimated mean of DLN
    var1.MOVER  = exp((2*mu_hat)+sigmaSq_hat)
    var2.MOVER  = delta_hat*(1-delta_hat)
    var3.MOVER = 0.5*(1-delta_hat)*((2*sigmaSq_hat)+sigmaSq_hat^2)
    var_thetahat.MOVER = var1.MOVER *(var2.MOVER +var3.MOVER)/n
    w_i = 1/var_thetahat.MOVER
    
    #The estimated mean of DLN
    Bessel1 = 1+((size_nonzero-1)/size_nonzero)*t
    Bessel2 = ((size_nonzero-1)^3/(factorial(2)*size_nonzero^2))*(t^2/(size_nonzero+1))
    Bessel3 = ((size_nonzero-1)^5/(factorial(3)*size_nonzero^3))*(t^3/((size_nonzero+1)*(size_nonzero+3)))
    Bessel4 = ((size_nonzero-1)^7/(factorial(4)*size_nonzero^4))*(t^4/((size_nonzero+1)*(size_nonzero+3)*(size_nonzero+5)))
    Bessel5 = ((size_nonzero-1)^9/(factorial(5)*size_nonzero^5))*(t^5/((size_nonzero+1)*(size_nonzero+3)*(size_nonzero+5)*(size_nonzero+7)))
    Bessel= Bessel1+Bessel2+Bessel3+Bessel4+Bessel5
    theta_i.MOVER = (1-delta_hat)*exp(mu_hat)*Bessel
    
    Temp.l.MOVER = (sum(w_i*theta_i.MOVER)/sum(w_i))- sqrt(sum(w_i^2*(theta_i.MOVER-l.MOVER)^2)/sum(w_i^2))      
    Temp.u.MOVER = (sum(w_i*theta_i.MOVER)/sum(w_i))+ sqrt(sum(w_i^2*(theta_i.MOVER-u.MOVER)^2)/sum(w_i^2))      
    
    l.mover = ifelse(w_i>0, Temp.l.MOVER, Temp.u.MOVER)
    u.mover = ifelse(w_i>0, Temp.u.MOVER, Temp.l.MOVER)
    
    L.mover = (sum(w_i*theta_i.MOVER)/sum(w_i))- sqrt(sum(w_i^2*(theta_i.MOVER-l.mover)^2)/sum(w_i^2))      
    U.mover = (sum(w_i*theta_i.MOVER)/sum(w_i))+ sqrt(sum(w_i^2*(theta_i.MOVER-u.mover)^2)/sum(w_i^2))     
    CP.mover[i] = ifelse(L.mover < theta && theta< U.mover,1,0)                  
    Length.mover[i] = U.mover - L.mover                                          
    
    
    #####################################################################
    #                  Algorithm 4: PB                                  #
    #####################################################################
    
    #Estimate SigmaSq
    w_hat.i = size_nonzero/sigmaSq_hat.MLE
    lnMeanDLN_hat = sum(w_hat.i*(mu_hat+(sigmaSq_hat.MLE/2)+log(1-delta_hat)))/sum(w_hat.i)
    AA = mu_hat - (lnMeanDLN_hat-log(1-delta_hat))
    
    #Estimate MLE of SigmaSq 
    sigmaSq_hat.MLEPB = 2*(sqrt(1 + sigmaSq_hat.MLE + AA^2)-1)
    w_hat.i.MLE = size_nonzero/sigmaSq_hat.MLEPB
    delta_hat.i.MLE = size_zero/(n+w_hat.i.MLE)
    lnMeanDLN_hat.MLE = sum(w_hat.i.MLE*(mu_hat+(sigmaSq_hat.MLEPB/2)+log(1-delta_hat.i.MLE)))/sum(w_hat.i.MLE)  
    
    h = 1
    while(h <= m){
      for(i_k in 1:k){
        x_DLN = x_data[which(x_data$xlabel==i_k),]
        
        A = sample(1:n[i_k],size=n[i_k],replace=TRUE)
        x_DLN.BS = x_DLN[A,]
        
        nonzero_value.BS = log(x_DLN.BS[which(x_DLN.BS$data!=0),'data'])
        size_nonzero.BS =length(nonzero_value.BS)
        size_zero.BS = n[i_k]-size_nonzero.BS
        mu_hat.BS = mean(nonzero_value.BS)
        sigmaSq_hat.BS = var(nonzero_value.BS)
        deltahat.BS = size_zero.BS/n[i_k]
        
        mu_hat.PB.i[i_k] = rnorm(1, mu_hat.BS, sqrt(sigmaSq_hat.BS/size_nonzero.BS))
        sigmaSq_hat.PB.i[i_k] = sigmaSq_hat.BS*rchisq(1,size_nonzero.BS-1)/(size_nonzero.BS-1)
        delta_hat.PB.i[i_k] = rbeta(1, size_zero.BS+0.5, size_nonzero.BS+0.5)
        size_nonzero.PB.i[i_k] = n[i_k]*(1-delta_hat.PB.i[i_k])
      }
      
      w.PB.i = size_nonzero.PB.i/sigmaSq_hat.PB.i
      lnMeanDLN.PB = sum(w.PB.i*(mu_hat.PB.i+(sigmaSq_hat.PB.i/2)+log(1-delta_hat.PB.i)))/sum(w.PB.i)
      
      AA.MLE = mu_hat.PB.i - (lnMeanDLN.PB-log(1-delta_hat.PB.i))
      sigmaSq_hat.MLE.PB = 2*(sqrt(1 + sigmaSq_hat.PB.i + AA.MLE^2)-1)
      w_hat.i.MLE.PB = size_nonzero.PB.i/sigmaSq_hat.MLE.PB
      delta_hat.i.MLE.PB = (n-size_nonzero.PB.i)/(n+w_hat.i.MLE.PB)
      lnMeanDLN_hat.MLE.PB = sum(w_hat.i.MLE.PB*(mu_hat.PB.i+(sigmaSq_hat.MLE.PB/2)+log(1-delta_hat.i.MLE.PB)))/sum(w_hat.i.MLE.PB)  
      
      Q.PB[h] = (lnMeanDLN_hat.MLE.PB-lnMeanDLN_hat.MLE)^2*sum(w_hat.i.MLE.PB)
      
      if(is.na(Q.PB[h])){next}
      h = h+1
    }
    
    q.PB = quantile(Q.PB, 1-alpha)
    L.PB = exp(lnMeanDLN_hat.MLE - sqrt(q.PB/sum(w_hat.i.MLE)))
    U.PB = exp(lnMeanDLN_hat.MLE + sqrt(q.PB/sum(w_hat.i.MLE)))
    CP.PB[i] = ifelse(L.PB < theta && theta < U.PB,1,0)
    Length.PB[i]= U.PB - L.PB
    
    
    #####################################################################
    #               Algorithm 5: HPD-JR and HPD-NGB                     #
    #####################################################################
    
    ############################### HPD-JR ##############################
    
    r = 0
    while(r < m){
      
      v = size_nonzero+1
      shap = v/2
      scal = v*sigmaSq_hat/2
      post.sigmaSq.JR = rinvgamma(k, shap, scal)
      
      Tem.post.mu.JR = numeric()
      for(k in 1:length(post.sigmaSq.JR)){
        Tem.post.mu.JR[k] = rnorm(1, mu_hat, sqrt(post.sigmaSq.JR[k]/size_nonzero))
      }
      post.mu.JR = Tem.post.mu.JR
      post.delta.JR  = rbeta(k, size_zero+0.5, size_nonzero+1.5)   
      postMeanDLN.JR = (1-post.delta.JR)*exp(post.mu.JR + post.sigmaSq.JR/2)
      
      #Variance of meanDLN_hat based on Owen and DeRouen 1980
      var1.JR = exp((2*post.mu.JR)+post.sigmaSq.JR)
      var2.JR = post.delta.JR*(1-post.delta.JR)
      var3.JR = 0.5*(1-post.delta.JR)*((2*post.sigmaSq.JR)+post.sigmaSq.JR^2)
      var_thetahat.JR = var1.JR*(var2.JR+var3.JR)/n
      
      w_i.JR = 1/var_thetahat.JR
      theta.JR [r] = sum(w_i.JR*postMeanDLN.JR)/sum(w_i.JR)
      
      r=r+1
    }
    
    HPD.JR = hdi(theta.JR, credMass = 1-alpha)
    L.HPD.JR  = quantile(HPD.JR ,0.025)                                       
    U.HPD.JR  = quantile(HPD.JR ,0.975)           
    CP.HPD.JR [i] = ifelse(L.HPD.JR < theta && theta < U.HPD.JR ,1,0)
    Length.HPD.JR [i] = U.HPD.JR -L.HPD.JR     
    
    
    ############################### HPD-NGB ##############################
    
    s = 0
    while(s < m){
      
      alpha_n = (size_nonzero-1)/2
      beta_n = (sum_nonzeroSq-(size_nonzero*mu_hat^2))/2
      post.sigmaSq.NGB = rinvgamma(k, alpha_n, beta_n)
      post.mu.NGB = rst(k, mu_hat, sqrt(sigmaSq_hat/size_nonzero), v)
      post.delta.NGB  = rbeta(k, size_zero+0.5, size_nonzero+0.5)   
      postMeanDLN.NGB  = (1-post.delta.NGB)*exp(post.mu.NGB + post.sigmaSq.NGB/2)
      
      #Variance of meanDLN_hat based on Owen and DeRouen 1980
      var1.NGB = exp((2*post.mu.NGB)+post.sigmaSq.NGB)
      var2.NGB = post.delta.NGB*(1-post.delta.NGB)
      var3.NGB = 0.5*(1-post.delta.NGB)*((2*post.sigmaSq.NGB)+post.sigmaSq.NGB^2)
      var_thetahat.NGB = var1.NGB*(var2.NGB+var3.NGB)/n
      
      w_i.NGB = 1/var_thetahat.NGB
      theta.NGB [s] = sum(w_i.NGB*postMeanDLN.NGB)/sum(w_i.NGB)
      s=s+1
    }
    
    HPD.NGB = hdi(theta.NGB, credMass = 1-alpha)
    L.HPD.NGB  = quantile(HPD.NGB ,0.025)                                       
    U.HPD.NGB  = quantile(HPD.NGB ,0.975)           
    CP.HPD.NGB [i] = ifelse(L.HPD.NGB < theta && theta < U.HPD.NGB ,1,0)
    Length.HPD.NGB [i] = U.HPD.NGB -L.HPD.NGB     
    
    
  }                                                                                #end loop i
  
  ACP.FGCI = mean(CP.FGCI)	                                                   
  Ex.length.FGCI = mean(Length.FGCI)                                           
  
  ACP.LS = mean(CP.LS)	                                                          
  Ex.length.LS = mean(Length.LS)                                                 
  
  ACP.mover = mean(CP.mover)	                                                
  Ex.length.mover = mean(Length.mover)                                       
  
  ACP.PB = mean(CP.PB)	                                                
  Ex.length.PB = mean(Length.PB)                                       
  
  
  ACP.HPD.JR = mean(CP.HPD.JR)	                                                
  Ex.length.HPD.JR = mean(Length.HPD.JR)                                       
  
  ACP.HPD.NGB = mean(CP.HPD.NGB)	                                                
  Ex.length.HPD.NGB = mean(Length.HPD.NGB)                                       
  
  
  elapseTime = proc.time()-starttime                                   
  
  cat("k =",k,"M =",M,"m =",m, "nTotal =",nTotal, "meanTotal=",meanTotal, "sigmaSqTotal=",sigmaSqTotal,"deltaTotal =",deltaTotal, "theta=",theta,"\n",
      
      "CP|FGCI. =",ACP.FGCI ,"Length|FGCI  =",Ex.length.FGCI,"\n",
      "CP|LS =",ACP.LS,"Length|LS =",Ex.length.LS,"\n",
      "CP|MOVER.W =",ACP.mover,"Length|MOVER.W =",Ex.length.mover,"\n",
      "CP|PB =",ACP.PB,"Length|PB =",Ex.length.PB,"\n",
      "CP|HPD.JR =",ACP.HPD.JR,"Length|HPD.JR =",Ex.length.HPD.JR,"\n",
      "CP|HPD.NGB =",ACP.HPD.NGB,"Length|HPD.NGB =",Ex.length.HPD.NGB)
  
  
  result = c("k =",k,"M =",M,"m =",m, "nTotal =",nTotal, "meanTotal=",meanTotal, "sigmaSqTotal=",sigmaSqTotal,"deltaTotal =",deltaTotal, "theta=",theta,
             
             
             "CP|FGCI =",ACP.FGCI,
             "CP|LS =",ACP.LS,
             "CP|MOVER.W =",ACP.mover,
             "CP|PB =",ACP.PB,
             "CP|HPD.JR =",ACP.HPD.JR,
             "CP|HPD.NGB =",ACP.HPD.NGB,
             
             
             "Length|FGCI  =",Ex.length.FGCI,
             "Length|LS =",Ex.length.LS,
             "Length|MOVER.W =",Ex.length.mover,
             "Length|PB =",Ex.length.PB,
             "Length|HPD.JR =",Ex.length.HPD.JR,
             "Length|HPD.NGB =",Ex.length.HPD.NGB,
             
             "Time =",elapseTime["elapsed"])
  return(result)
  
  
}





##################################################################################
# k = 10 sample cases                                                             #
##################################################################################

xGen.k10 = function(k,nTotal,Ex_LN,cv_LN,deltaTotal){
  sum_nonzeroSq = rep(0,k)
  mu_hat = rep(0,k)
  sigmaSq_hat= rep(0,k)
  sigmaSq_hat.MLE= rep(0,k)
  delta_hatnon= rep(0,k)
  delta_hat= rep(0,k)
  size_non = rep(0,k)
  size_zero = rep(0,k)
  xlist = list()
  x_data = data.frame()
  
  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)])                         #find values is not zero
      Temp.size_nonzero =length(Temp.nonzero_value)                          #find size of nonzero
      if(Temp.size_nonzero < (nTotal[l]-nTotal[l]*deltaTotal[l])){
        next
      }
      nonzero_value = Temp.nonzero_value
      size_nonzero = Temp.size_nonzero
      sum_nonzeroSq[l] = sum(nonzero_value^2)
      mu_hat[l] = mean(nonzero_value)
      sigmaSq_hat[l] = var(nonzero_value)
      
      sigmaSq_hat.MLE[l] = (1/size_nonzero)*(sum(nonzero_value^2)-(size_nonzero*mu_hat[l]^2))
      
      delta_hatnon[l] = size_nonzero/nTotal[l]
      delta_hat[l] = 1- delta_hatnon[l]
      size_non[l] = size_nonzero
      size_zero[l] = nTotal[l]-size_nonzero
      
      tmp_x = data.frame(x)
      tmp_x$xlabel = l
      names(tmp_x) = c("data","xlabel")
      x_data = rbind(x_data,tmp_x)
      
      break
    }
  }
  codeN = paste(nTotal, collapse = '')
  filename = paste("result_x_K",k,"_",codeN,".csv", sep='')
  write.csv(x_data,filename,row.names = FALSE)
  return(cbind(size_non,size_zero, mu_hat, sigmaSq_hat, sigmaSq_hat.MLE, delta_hat, sum_nonzeroSq))	
}

CI.CommonMean.k10 = function(M=5000,m=2500,k=10, 
                             n1=30, n2=30, n3=30, n4=30, n5=30, n6=50, n7=50, n8=50, n9=50, n10=50,
                             
                             delta1=0.1, delta2=0.1, delta3=0.1, delta4=0.1, delta5=0.1, 
                             delta6=0.2, delta7=0.2, delta8=0.2, delta9=0.2, delta10=0.2,
                             
                             sigmaSq1=1, sigmaSq2=1, sigmaSq3=1, sigmaSq4=1, sigmaSq5=1, 
                             sigmaSq6=2, sigmaSq7=2, sigmaSq8=2, sigmaSq9=2, sigmaSq10=2,
                             alpha=0.05){
  starttime = proc.time()
  
  
  nTotal = c(n1,n2,n3,n4,n5,n6,n7,n8,n9,n10)
  sigmaSqTotal = c(sigmaSq1,sigmaSq2,sigmaSq3,sigmaSq4,sigmaSq5,sigmaSq6,sigmaSq7,sigmaSq8,sigmaSq9,sigmaSq10)
  meanTotal =  rep(0, times = 10)
  deltaTotal = c(delta1,delta2,delta3,delta4,delta5,delta6,delta7,delta8,delta9,delta10)
  size_nonzeroTotal = nTotal*(1-deltaTotal)
  k = length(nTotal)
  
  Ex_LN = exp(meanTotal+(sigmaSqTotal/2))
  cv_LN = sqrt(exp(sigmaSqTotal)-1)
  
  var1 = exp((2*meanTotal)+sigmaSqTotal)
  var2 = deltaTotal*(1-deltaTotal)
  var3 = 0.5*(1-deltaTotal)*((2*sigmaSqTotal)+sigmaSqTotal^2)
  Var_theta = var1*(var2+var3)/nTotal
  
  w_i = 1/Var_theta
  theta = sum(w_i*(1-deltaTotal)*exp(meanTotal+(sigmaSqTotal/2)))/sum(w_i)
  
  #FGCI
  theta.FGCI = matrix(rep(0,m))
  L.FGCI = numeric()
  U.FGCI = numeric()
  CP.FGCI = rep(0,M)
  Length.FGCI = rep(0,M)
  
  #LS
  L.LS = numeric()
  U.LS  = numeric()
  CP.LS  = rep(0,M) 
  Length.LS  = rep(0,M)
  
  #MOVER
  L.mover = numeric()
  U.mover = numeric()
  CP.mover = rep(0,M)
  Length.mover = rep(0,M)
  
  #PB
  mu_hat.PB.i = rep(0,k)
  sigmaSq_hat.PB.i = rep(0,k)
  delta_hat.PB.i = rep(0,k)
  size_nonzero.PB.i = rep(0,k)
  Q.PB = rep(0,m)
  L.PB = numeric()
  U.PB = numeric()
  CP.PB = rep(0,M)
  Length.PB = rep(0,M)
  
  #HPD-JR
  post.mu.JR = rep(0,k)
  post.sigmaSq.JR = rep(0,k)
  post.delta.JR = rep(0,k)
  postMeanDLN.JR = rep(0,k)
  theta.JR = rep(0,m)
  L.HPD.JR = numeric()
  U.HPD.JR = numeric()
  CP.HPD.JR = rep(0,M)
  Length.HPD.JR = rep(0,M)
  
  #HPD-NGB
  post.mu.NGB = rep(0,k)
  post.sigmaSq.NGB = rep(0,k)
  post.delta.NGB = rep(0,k)
  postMeanDLN.NGB = rep(0,k) 
  theta.NGB = rep(0,m)
  L.HPD.NGB = numeric()
  U.HPD.NGB = numeric()
  CP.HPD.NGB = rep(0,M)
  Length.HPD.NGB = rep(0,M)
  
  
  i = 0
  while(i < M){                                                   # begin loop i
    
    
    data = xGen.k10(k, nTotal, Ex_LN, cv_LN,deltaTotal)
    codeN = paste(nTotal, collapse = '')
    filename = paste("result_x_K",k,"_",codeN,".csv", sep='')
    x_data = read.csv(filename)
    
    
    size_nonzero = data[,1]
    size_zero = data[,2]
    mu_hat = data[,3]
    sigmaSq_hat  = data[,4]
    sigmaSq_hat.MLE  = data[,5]
    delta_hat= data[,6]
    sum_nonzeroSq= data[,7]
    n = nTotal
    
    if(sum(delta_hat==0)>0){next}
    i=i+1
    if(i%%1000==0){
      print(i)
    }
    
    
    ##############################################################################
    #                         Algorithm 1: FGCI                                  #
    ##############################################################################
    e = 1
    while(e <= m){
      
      USq = rchisq(k,size_nonzero-1)/(size_nonzero-1)
      mu_FGPQ = mu_hat - ((rnorm(k)/sqrt(USq))*sqrt(sigmaSq_hat/size_nonzero))
      sigmaSq_FGPQ = sigmaSq_hat/USq
      deltanon_FGPQ = rbeta(k,size_nonzero+0.5, size_zero+0.5)
      delta_FGPQ = 1-deltanon_FGPQ
      MeanDLN_FGPQ_i  = deltanon_FGPQ*exp(mu_FGPQ+(sigmaSq_FGPQ/2))
      
      #Variance of meanDLN_hat based on Owen and DeRouen 1980
      var1.FGCI = exp((2*mu_FGPQ)+sigmaSq_FGPQ)
      var2.FGCI = delta_FGPQ*(1-delta_FGPQ)
      var3.FGCI = 0.5*(1-delta_FGPQ)*((2*sigmaSq_FGPQ)+sigmaSq_FGPQ^2)
      Var_thetahat.FGCI = var1.FGCI*(var2.FGCI+var3.FGCI)/n
      
      w_i.FGCI = 1/Var_thetahat.FGCI
      theta.FGCI[e] = sum(w_i.FGCI*MeanDLN_FGPQ_i)/sum(w_i.FGCI)
      e=e+1
    }
    
    L.FGCI = quantile(theta.FGCI,0.025)                                       
    U.FGCI = quantile(theta.FGCI,0.975)      
    CP.FGCI[i] = ifelse(L.FGCI< theta && theta <U.FGCI,1,0)
    Length.FGCI[i] = U.FGCI - L.FGCI
    
    #####################################################################
    #                     Algorithm 2: LS                               #
    #####################################################################
    
    k_df =size_nonzero-1
    
    #Bessel function
    t = sigmaSq_hat/2
    
    Bessel1 = 1+((size_nonzero-1)/size_nonzero)*t
    Bessel2 = ((size_nonzero-1)^3/(factorial(2)*size_nonzero^2))*(t^2/(size_nonzero+1))
    Bessel3 = ((size_nonzero-1)^5/(factorial(3)*size_nonzero^3))*(t^3/((size_nonzero+1)*(size_nonzero+3)))
    Bessel4 = ((size_nonzero-1)^7/(factorial(4)*size_nonzero^4))*(t^4/((size_nonzero+1)*(size_nonzero+3)*(size_nonzero+5)))
    Bessel5 = ((size_nonzero-1)^9/(factorial(5)*size_nonzero^5))*(t^5/((size_nonzero+1)*(size_nonzero+3)*(size_nonzero+5)*(size_nonzero+7)))
    Bessel= Bessel1+Bessel2+Bessel3+Bessel4+Bessel5
    theta_hat_i = (1-delta_hat)*exp(mu_hat)*Bessel
    
    #Variance of meanDLN_hat based on Owen and DeRouen 1980
    var1.LS  = exp((2*mu_hat)+sigmaSq_hat)
    var2.LS  = delta_hat*(1-delta_hat)
    var3.LS = 0.5*(1-delta_hat)*((2*sigmaSq_hat)+sigmaSq_hat^2)
    var_thetahat.LS = var1.LS *(var2.LS +var3.LS)/n
    
    w_i = 1/var_thetahat.LS
    theta_hat = sum(w_i*theta_hat_i)/sum(w_i)
    z =  qnorm(1-(alpha/2),0,1)
    sd.LS = sqrt(1/(sum(w_i)))
    
    L.LS = theta_hat-(z*sd.LS)                                      
    U.LS = theta_hat+(z*sd.LS)           
    CP.LS[i] = ifelse(L.LS< theta && theta<U.LS,1,0)                 
    Length.LS[i] = U.LS-L.LS                                        
    
    ##############################################################################
    #                          Algorithm 3:  MOVER                               #
    ##############################################################################
    
    #Wilson interval for delta
    r = qnorm(1-(alpha/2))
    deltanon_hat = 1 - delta_hat
    CI.delta1 = deltanon_hat +(0.5*r^2/n)
    CI.delta2 = sqrt((deltanon_hat*(1-deltanon_hat)+ (0.25*r^2/n))/n)
    CI.delta3 = 1+(r^2/n) 
    l.delta = log((CI.delta1-(r*CI.delta2))/CI.delta3)
    u.delta = log((CI.delta1+(r*CI.delta2))/CI.delta3)
    
    
    #CI for sigmaSq
    df = size_nonzero-1
    locate_L = qchisq(alpha/2,df) 
    locate_U = qchisq(1-(alpha/2),df) 
    
    #MOVER CI for mu + sigmaSq/2
    sd_mu = qnorm(1-(alpha/2))*sqrt(sigmaSq_hat/size_nonzero)
    sd_sigmaSq_L = (sigmaSq_hat/2)*(1-((size_nonzero-1)/(locate_U)))
    sd_sigmaSq_U = (sigmaSq_hat/2)*(((size_nonzero-1)/(locate_L))-1)
    L_the2_3= (mu_hat+sigmaSq_hat/2)-sqrt(sd_mu^2 + sd_sigmaSq_L^2)
    U_the2_3= (mu_hat+sigmaSq_hat/2)+sqrt(sd_mu^2 + sd_sigmaSq_U^2)
    
    #MOVER CI for log(1-delta_hat) + mu_hat + sigmaSq_hat/2
    l.MOVER = exp((log(1-delta_hat)+(mu_hat + sigmaSq_hat/2))-sqrt((log(1-delta_hat)-l.delta)^2+((mu_hat + sigmaSq_hat/2)-L_the2_3)^2))
    u.MOVER = exp((log(1-delta_hat)+(mu_hat + sigmaSq_hat/2))+sqrt((u.delta-log(1-delta_hat))^2+(U_the2_3-(mu_hat + sigmaSq_hat/2))^2))
    
    #Variance of estimated mean of DLN
    var1.MOVER  = exp((2*mu_hat)+sigmaSq_hat)
    var2.MOVER  = delta_hat*(1-delta_hat)
    var3.MOVER = 0.5*(1-delta_hat)*((2*sigmaSq_hat)+sigmaSq_hat^2)
    var_thetahat.MOVER = var1.MOVER *(var2.MOVER +var3.MOVER)/n
    w_i = 1/var_thetahat.MOVER
    
    #The estimated mean of DLN
    t = sigmaSq_hat/2
    Bessel1 = 1+((size_nonzero-1)/size_nonzero)*t
    Bessel2 = ((size_nonzero-1)^3/(factorial(2)*size_nonzero^2))*(t^2/(size_nonzero+1))
    Bessel3 = ((size_nonzero-1)^5/(factorial(3)*size_nonzero^3))*(t^3/((size_nonzero+1)*(size_nonzero+3)))
    Bessel4 = ((size_nonzero-1)^7/(factorial(4)*size_nonzero^4))*(t^4/((size_nonzero+1)*(size_nonzero+3)*(size_nonzero+5)))
    Bessel5 = ((size_nonzero-1)^9/(factorial(5)*size_nonzero^5))*(t^5/((size_nonzero+1)*(size_nonzero+3)*(size_nonzero+5)*(size_nonzero+7)))
    Bessel= Bessel1+Bessel2+Bessel3+Bessel4+Bessel5
    theta_i.MOVER = (1-delta_hat)*exp(mu_hat)*Bessel
    
    Temp.l.MOVER = (sum(w_i*theta_i.MOVER)/sum(w_i))- sqrt(sum(w_i^2*(theta_i.MOVER-l.MOVER)^2)/sum(w_i^2))      
    Temp.u.MOVER = (sum(w_i*theta_i.MOVER)/sum(w_i))+ sqrt(sum(w_i^2*(theta_i.MOVER-u.MOVER)^2)/sum(w_i^2))     
    
    l.mover = ifelse(w_i>0, Temp.l.MOVER, Temp.u.MOVER)
    u.mover = ifelse(w_i>0, Temp.u.MOVER, Temp.l.MOVER)
    
    L.mover = (sum(w_i*theta_i.MOVER)/sum(w_i))- sqrt(sum(w_i^2*(theta_i.MOVER-l.mover)^2)/sum(w_i^2))      
    U.mover = (sum(w_i*theta_i.MOVER)/sum(w_i))+ sqrt(sum(w_i^2*(theta_i.MOVER-u.mover)^2)/sum(w_i^2))      
    CP.mover[i] = ifelse(L.mover < theta && theta< U.mover,1,0)                 
    Length.mover[i] = U.mover - L.mover                                          
    
    
    #####################################################################
    #                  Algorithm 4: PB                                  #
    #####################################################################
    #Estimate SigmaSq
    w_hat.i = size_nonzero/sigmaSq_hat.MLE
    lnMeanDLN_hat = sum(w_hat.i*(mu_hat+(sigmaSq_hat.MLE/2)+log(1-delta_hat)))/sum(w_hat.i)
    AA = mu_hat - (lnMeanDLN_hat-log(1-delta_hat))
    
    #Estimate MLE of SigmaSq
    sigmaSq_hat.MLEPB = 2*(sqrt(1 + sigmaSq_hat.MLE + AA^2)-1)
    w_hat.i.MLE = size_nonzero/sigmaSq_hat.MLEPB
    delta_hat.i.MLE = size_zero/(n+w_hat.i.MLE)
    lnMeanDLN_hat.MLE = sum(w_hat.i.MLE*(mu_hat+(sigmaSq_hat.MLEPB/2)+log(1-delta_hat.i.MLE)))/sum(w_hat.i.MLE)  
    
    h = 1
    while(h <= m){
      for(i_k in 1:k){
        x_DLN = x_data[which(x_data$xlabel==i_k),]
        
        A = sample(1:n[i_k],size=n[i_k],replace=TRUE)
        x_DLN.BS = x_DLN[A,]
        
        nonzero_value.BS = log(x_DLN.BS[which(x_DLN.BS$data!=0),'data'])
        size_nonzero.BS =length(nonzero_value.BS)
        size_zero.BS = n[i_k]-size_nonzero.BS
        mu_hat.BS = mean(nonzero_value.BS)
        sigmaSq_hat.BS = var(nonzero_value.BS)
        deltahat.BS = size_zero.BS/n[i_k]
        
        mu_hat.PB.i[i_k] = rnorm(1, mu_hat.BS, sqrt(sigmaSq_hat.BS/size_nonzero.BS))
        sigmaSq_hat.PB.i[i_k] = sigmaSq_hat.BS*rchisq(1,size_nonzero.BS-1)/(size_nonzero.BS-1)
        delta_hat.PB.i[i_k] = rbeta(1, size_zero.BS+0.5, size_nonzero.BS+0.5)
        size_nonzero.PB.i[i_k] = n[i_k]*(1-delta_hat.PB.i[i_k])
        
        
      }
      
      w.PB.i = size_nonzero.PB.i/sigmaSq_hat.PB.i
      lnMeanDLN.PB = sum(w.PB.i*(mu_hat.PB.i+(sigmaSq_hat.PB.i/2)+log(1-delta_hat.PB.i)))/sum(w.PB.i)
      
      AA.MLE = mu_hat.PB.i - (lnMeanDLN.PB-log(1-delta_hat.PB.i))
      sigmaSq_hat.MLE.PB = 2*(sqrt(1 + sigmaSq_hat.PB.i + AA.MLE^2)-1)
      w_hat.i.MLE.PB = size_nonzero.PB.i/sigmaSq_hat.MLE.PB
      delta_hat.i.MLE.PB = (n-size_nonzero.PB.i)/(n+w_hat.i.MLE.PB)
      lnMeanDLN_hat.MLE.PB = sum(w_hat.i.MLE.PB*(mu_hat.PB.i+(sigmaSq_hat.MLE.PB/2)+log(1-delta_hat.i.MLE.PB)))/sum(w_hat.i.MLE.PB)  
      
      Q.PB[h] = (lnMeanDLN_hat.MLE.PB-lnMeanDLN_hat.MLE)^2*sum(w_hat.i.MLE.PB)
      
      if(is.na(Q.PB[h])){next}
      h = h+1
    }
    
    q.PB = quantile(Q.PB, 1-alpha)
    L.PB = exp(lnMeanDLN_hat.MLE - sqrt(q.PB/sum(w_hat.i.MLE)))
    U.PB = exp(lnMeanDLN_hat.MLE + sqrt(q.PB/sum(w_hat.i.MLE)))
    CP.PB[i] = ifelse(L.PB < theta && theta < U.PB,1,0)
    Length.PB[i]= U.PB - L.PB
    
    
    #####################################################################
    #               Algorithm 5: HPD-JR and HPD-NGB                     #
    #####################################################################
    
    ############################### HPD-JR ##############################
    
    r = 0
    while(r < m){
      
      v = size_nonzero+1
      shap = v/2
      scal = v*sigmaSq_hat/2
      post.sigmaSq.JR = rinvgamma(k, shap, scal)
      
      Tem.post.mu.JR = numeric()
      for(k in 1:length(post.sigmaSq.JR)){
        Tem.post.mu.JR[k] = rnorm(1, mu_hat, sqrt(post.sigmaSq.JR[k]/size_nonzero))
      }
      post.mu.JR = Tem.post.mu.JR
      post.delta.JR  = rbeta(k, size_zero+0.5, size_nonzero+1.5)   
      postMeanDLN.JR = (1-post.delta.JR)*exp(post.mu.JR + post.sigmaSq.JR/2)
      
      #Variance of meanDLN_hat based on Owen and DeRouen 1980
      var1.JR = exp((2*post.mu.JR)+post.sigmaSq.JR)
      var2.JR = post.delta.JR*(1-post.delta.JR)
      var3.JR = 0.5*(1-post.delta.JR)*((2*post.sigmaSq.JR)+post.sigmaSq.JR^2)
      var_thetahat.JR = var1.JR*(var2.JR+var3.JR)/n
      
      w_i.JR = 1/var_thetahat.JR
      theta.JR [r] = sum(w_i.JR*postMeanDLN.JR)/sum(w_i.JR)
      
      r=r+1
    }
    
    HPD.JR = hdi(theta.JR, credMass = 1-alpha)
    L.HPD.JR  = quantile(HPD.JR ,0.025)                                       
    U.HPD.JR  = quantile(HPD.JR ,0.975)           
    CP.HPD.JR [i] = ifelse(L.HPD.JR < theta && theta < U.HPD.JR ,1,0)
    Length.HPD.JR [i] = U.HPD.JR -L.HPD.JR     
    
    
    ############################### HPD-NGB ##############################
    
    s = 0
    while(s < m){
      
      alpha_n = (size_nonzero-1)/2
      beta_n = (sum_nonzeroSq-(size_nonzero*mu_hat^2))/2
      post.sigmaSq.NGB = rinvgamma(k, alpha_n, beta_n)
      post.mu.NGB = rst(k, mu_hat, sqrt(sigmaSq_hat/size_nonzero), v)
      post.delta.NGB  = rbeta(k, size_zero+0.5, size_nonzero+0.5)   
      postMeanDLN.NGB  = (1-post.delta.NGB)*exp(post.mu.NGB + post.sigmaSq.NGB/2)
      
      #Variance of meanDLN_hat based on Owen and DeRouen 1980
      var1.NGB = exp((2*post.mu.NGB)+post.sigmaSq.NGB)
      var2.NGB = post.delta.NGB*(1-post.delta.NGB)
      var3.NGB = 0.5*(1-post.delta.NGB)*((2*post.sigmaSq.NGB)+post.sigmaSq.NGB^2)
      var_thetahat.NGB = var1.NGB*(var2.NGB+var3.NGB)/n
      
      w_i.NGB = 1/var_thetahat.NGB
      theta.NGB [s] = sum(w_i.NGB*postMeanDLN.NGB)/sum(w_i.NGB)
      
      s=s+1
    }
    
    HPD.NGB = hdi(theta.NGB, credMass = 1-alpha)
    L.HPD.NGB  = quantile(HPD.NGB ,0.025)                                       
    U.HPD.NGB  = quantile(HPD.NGB ,0.975)           
    CP.HPD.NGB [i] = ifelse(L.HPD.NGB < theta && theta < U.HPD.NGB ,1,0)
    Length.HPD.NGB [i] = U.HPD.NGB -L.HPD.NGB     
    
    
  }                                                                                #end loop i
  
  ACP.FGCI = mean(CP.FGCI)	                                                    
  Ex.length.FGCI = mean(Length.FGCI)                                          
  
  ACP.LS = mean(CP.LS)	                                                          
  Ex.length.LS = mean(Length.LS)                                                
  
  ACP.mover = mean(CP.mover)	                                                
  Ex.length.mover = mean(Length.mover)                                       
  
  ACP.PB = mean(CP.PB)	                                                
  Ex.length.PB = mean(Length.PB)                                      
  
  ACP.HPD.JR = mean(CP.HPD.JR)	                                                
  Ex.length.HPD.JR = mean(Length.HPD.JR)                                       
  
  ACP.HPD.NGB = mean(CP.HPD.NGB)	                                                
  Ex.length.HPD.NGB = mean(Length.HPD.NGB)                                       
  
  
  elapseTime = proc.time()-starttime                                   
  
  cat("k =",k,"M =",M,"m =",m, "nTotal =",nTotal, "meanTotal=",meanTotal, "sigmaSqTotal=",sigmaSqTotal,"deltaTotal =",deltaTotal, "theta=",theta,"\n",
      
      "CP|FGCI. =",ACP.FGCI ,"Length|FGCI  =",Ex.length.FGCI,"\n",
      "CP|LS =",ACP.LS,"Length|LS =",Ex.length.LS,"\n",
      "CP|mover =",ACP.mover,"Length|mover =",Ex.length.mover,"\n",
      "CP|PB =",ACP.PB,"Length|PB =",Ex.length.PB,"\n",
      "CP|HPD.JR =",ACP.HPD.JR,"Length|HPD.JR =",Ex.length.HPD.JR,"\n",
      "CP|HPD.NGB =",ACP.HPD.NGB,"Length|HPD.NGB =",Ex.length.HPD.NGB)
  
  
  result = c("k =",k,"M =",M,"m =",m, "nTotal =",nTotal, "meanTotal=",meanTotal, "sigmaSqTotal=",sigmaSqTotal,"deltaTotal =",deltaTotal, "theta=",theta,
             
             
             "CP|FGCI =",ACP.FGCI,
             "CP|LS =",ACP.LS,
             "CP|mover =",ACP.mover,
             "CP|PB =",ACP.PB,
             "CP|HPD.JR =",ACP.HPD.JR,
             "CP|HPD.NGB =",ACP.HPD.NGB,
             
             "Length|FGCI  =",Ex.length.FGCI,
             "Length|LS =",Ex.length.LS,
             "Length|mover =",Ex.length.mover,
             "Length|PB =",Ex.length.PB,
             "Length|HPD.JR =",Ex.length.HPD.JR,
             "Length|HPD.NGB =",Ex.length.HPD.NGB,
             
             "Time =",elapseTime["elapsed"])
  return(result)
  
  
}



##############################################################################
# Dataset 1:                                                                 #
##############################################################################
# R.code to apply with rainfall data on Augest 5, 2019                       #
##############################################################################
rm(list=ls())
library(MASS)

#Daily rainfall data in five Thailand’s regions 
#Northern
N = c(3.0,  2.6,  1.0,  3.6,  0.0, 13.2, 22.4,  1.4, 18.3,  0.0,
      15.5,  0.0,  0.0,  0.0,  9.8, 24.3, 24.6,  8.8,  0.0, 19.8,
      5.0, 12.3,  8.1,  4.8,  5.8, 17.0, 25.1,  8.3, 22.9, 26.9,
      0.0,  0.0,  5.0, 23.8, 16.0, 11.5,  1.2, 10.3,  1.7,  5.5,
      7.3, 24.3, 27.2, 12.6, 22.7, 0.0, 2.6,  0.0,  3.2,  2.6,
      2.0,  8.0,  1.9,  0.8,  2.2,  6.5,  0.0,  2.2, 0.0,  4.3,  0.2,  0.0)

#Northeastern
NE = c(3.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  1.7,  2.3,  0.5,
       3.9,  6.9,  2.2,  3.2,  5.3, 11.0,  0.6,  0.0,  1.0,  2.4, 13.2,  0.4,  0.0,
       1.3, 10.0,  2.5,  4.6,  0.0,  0.0, 40.0,  3.5,  0.0, 12.0, 15.0,  0.0,  1.5,
       0.7,  0.0,  3.0,  0.0,  0.0,  0.0, 29.4, 48.0,  0.0, 70.8,  3.5, 14.2,  7.0,
       0.0,  0.0,  0.0,  0.0, 10.8,  0.0,  6.3,  0.0,  4.0, 19.3,  0.0,  1.5, 18.5,
       42.0,  9.1,  6.0,  7.5,  0.0,  6.3,  0.0,  0.4,  0.0,  0.0,  0.0,  1.8,  0.0,
       0.0, 14.3,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0, 10.1,  0.0,  4.8,
       0.0,  0.0, 49.5, 10.5, 60.4, 12.7,  6.8, 69.3, 36.5,  8.6,  0.0,  0.0,  0.0,
       3.8,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0, 21.5,  2.5,  0.0,
       13.0, 26.2,  2.2,  3.0, 10.5,  0.0,  0.0,  0.0,  0.0,  4.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,  6.0,  0.0,  0.0,
       0.0,  2.3,  0.0,  1.0,  0.0,  0.0,  0.0,  4.6,  0.0, 10.0,  0.0,  9.5,  0.0,
       0.0,  0.0,  0.0, 20.3,  0.0,  2.4,  0.0,  0.0,  0.0,  0.0,  0.0,  3.2,  0.0,
       0.0,  0.0,  0.0,  0.0,  0.0,  4.8,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  5.4,
       0.0,  0.0, 12.0,  0.0,  0.0,  0.0, 11.0,  0.0,  0.0,  0.0,  0.3,  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)

#Central 
C = c(2.9,  0.2,  0.3,  2.5,  0.4,  0.4,  1.1,  0.0,  1.3,  0.1,  2.9,  0.0,  1.0,
      4.7,  0.5,  5.0,  2.5,  0.0,  0.0,  0.0,  0.0,  6.6,  0.0,  9.5,  5.1, 12.5,
      0.0,  0.0,  3.2,  0.0,  2.2,  3.2,  0.0,  0.0,  4.7, 19.3,  3.1,  2.9,  5.7,
      0.9,  0.0,  0.0,  2.6, 17.0,  0.0,  3.5,  0.0,  0.0,  0.0,  5.1, 60.4,  6.9,
      3.0, 15.1,  6.0, 13.4,  6.2)

#Eastern
E = c( 0.0,   3.2,  10.4,   1.1,   0.2,   4.3,   0.0,   0.0,   0.0,   0.0,   0.2,
       0.1,  62.8,  36.7,  15.6,  50.0,  35.5,  35.0,   5.9,   0.0,   0.0,   3.0,
       60.4,  60.0,  76.0,  79.7,  65.7, 108.0,  10.5)

#Southern
S = c(4.1,  0.0, 11.5,  2.5,  9.7, 10.4,  9.6, 19.0,  8.3,  0.0,  0.0,  0.0,  0.0, 17.8, 12.3,  
      2.5,  0.0,  0.9,  0.0,  2.6,  0.0,  0.0,  2.0,  0.0,  0.0,  0.0,  3.5,  0.0,  0.0,
      0.0,  0.0,  0.0,  0.0,  3.5, 13.6,  0.0,  0.0,  0.0,  0.0,  0.0,  4.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,  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.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.0,  0.0,
      0.0,  0.0, 36.1, 41.8, 30.0,  0.0,  2.7,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  
      6.2,  0.0,  0.0,  6.1,  0.0,  0.0,  0.0,  0.3,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,
      0.0,  0.0,  0.0)


##############################################################################
# R.code to construct histogram plots                                        #
##############################################################################
#Figure 4

png("Figure4.png", width = 6, height = 6, units = 'in', res = 300)
layout(matrix(c(1,1,2,2,3,3,0,4,4,5,5,0), 2, 6, byrow = TRUE))
hist(N, main="A", xlab = "Daily rainfall amounts",ylab = "Records",border = TRUE)
hist(NE, main="B", xlab = "Daily rainfall amounts",ylab = "Records")
hist(C, main="C", xlab = "Daily rainfall amounts",ylab = "Records")
hist(E, main="D", xlab = "Daily rainfall amounts",ylab = "Records")
hist(S, main="E", xlab = "Daily rainfall amounts",ylab = "Records")
dev.off() 
 

##############################################################################
# R.code to construct Normal Q-Q plots                                        #
##############################################################################
#Figure 6

png("Figure5.png", width = 6, height = 6, units = 'in', res = 300)
layout(matrix(c(1,1,2,2,3,3,0,4,4,5,5,0), 2, 6, byrow = TRUE))
log_N = log(N[N !=0])
log_NE = log(NE[NE !=0])
log_C = log(C[C !=0])
log_E = log(E[E !=0])
log_S = log(S[S !=0])

qqnorm(log_N, main = "A")
qqline(log_N,distribution = qnorm,probs = c(0.25, 0.75), qtype = 7)

qqnorm(log_NE, main = "B")
qqline(log_NE,distribution = qnorm, probs = c(0.25, 0.75), qtype = 7)

qqnorm(log_C, main = "C")
qqline(log_C,distribution = qnorm, probs = c(0.1, 0.85), qtype = 7)

qqnorm(log_E, main = "D")
qqline(log_E,distribution = qnorm, probs = c(0.20, 0.75), qtype = 7)

qqnorm(log_S, main = "E")
qqline(log_S,distribution = qnorm, probs = c(0.25, 0.75), qtype = 7)
dev.off()

##############################################################################
# R.code to estimate AIC results                                             #
##############################################################################

#Dataset: north
North_nonzero = N[N!=0]
cau_N = fitdistr(North_nonzero, "cauchy")
lognormal_N = fitdistr(North_nonzero, "lognormal")
normal_N = fitdistr(North_nonzero, "normal")
logistic_N = fitdistr(North_nonzero, "logistic")
t_N = fitdistr(North_nonzero, "t",df=length(North_nonzero)-1)

AIC(cau_N)
AIC(logistic_N)
AIC(lognormal_N)
AIC(normal_N)
AIC(t_N)

#Dataset: northeast
NE_nonzero = NE[NE !=0]
cau_NE = fitdistr(NE_nonzero, "cauchy")
lognormal_NE = fitdistr(NE_nonzero, "lognormal")
normal_NE = fitdistr(NE_nonzero, "normal")
logistic_NE = fitdistr(NE_nonzero, "logistic")
t_NE = fitdistr(NE_nonzero, "t",df=length(NE_nonzero)-1)

AIC(cau_NE)
AIC(logistic_NE)
AIC(lognormal_NE)
AIC(normal_NE)
AIC(t_NE)

#Dataset: central
C_nonzero = C[C !=0]
cau_C = fitdistr(C_nonzero, "cauchy")
lognormal_C = fitdistr(C_nonzero, "lognormal")
normal_C = fitdistr(C_nonzero, "normal")
logistic_C = fitdistr(C_nonzero, "logistic")
t_C = fitdistr(C_nonzero, "t",df=length(C_nonzero)-1)

AIC(cau_C)
AIC(logistic_C)
AIC(lognormal_C)
AIC(normal_C)
AIC(t_C)

#Dataset: east
E_nonzero = E[E!=0]
cau_E = fitdistr(E_nonzero, "cauchy")
lognormal_E = fitdistr(E_nonzero, "lognormal")
normal_E = fitdistr(E_nonzero, "normal")
logistic_E = fitdistr(E_nonzero, "logistic")
t_E = fitdistr(E_nonzero, "t", df=length(E_nonzero)-1)

AIC(cau_E)
AIC(logistic_E)
AIC(lognormal_E)
AIC(normal_E)
AIC(t_E)

#Dataset: south
S_nonzero = S[S!=0]
cau_S = fitdistr(S_nonzero, "cauchy")
lognormal_S = fitdistr(S_nonzero, "lognormal")
normal_S = fitdistr(S_nonzero, "normal")
logistic_S = fitdistr(S_nonzero, "logistic")
t_S = fitdistr(S_nonzero, "t",df=length(S_nonzero)-1)

AIC(cau_S)
AIC(logistic_S)
AIC(lognormal_S)
AIC(normal_S)
AIC(t_S)

AIC.N = data.frame()
AIC.NE = data.frame()
AIC.C = data.frame()
AIC.E = data.frame()
AIC.S = data.frame()
AIC.N = c(AIC(cau_N), AIC(logistic_N), AIC(lognormal_N), AIC(normal_N), AIC(t_N))
AIC.NE = c(AIC(cau_NE), AIC(logistic_NE), AIC(lognormal_NE), AIC(normal_NE), AIC(t_NE))
AIC.C = c(AIC(cau_C), AIC(logistic_C), AIC(lognormal_C), AIC(normal_C), AIC(t_C))
AIC.E = c(AIC(cau_E), AIC(logistic_E), AIC(lognormal_E), AIC(normal_E), AIC(t_E))
AIC.S = c(AIC(cau_S), AIC(logistic_S), AIC(lognormal_S), AIC(normal_S), AIC(t_S))
AIC = rbind(AIC.N, AIC.NE, AIC.C, AIC.E, AIC.S)

colnames(AIC) = c("Cauchy", "Logistic", "Lognormal", "Normal", "t")
AIC

##############################################################################
# Dataset 2:                                                                 #
##############################################################################
# R.code to apply with rainfall data on Augest 9, 2019                       #
##############################################################################
rm(list=ls())

#Daily rainfall data in five Thailand’s regions on Augest 9, 2019 
#Northern
N = c(9.5,  4.9,  0.0,  4.7,  0.0, 63.2,  9.6, 10.7, 13.0,  0.0,  0.0,  0.0,  0.0,  0.0,  2.8, 11.3,  0.6, 36.1,  
      0.0,  2.6,  5.0, 13.4, 12.3, 25.8,  30.2, 16.4,  6.0, 33.1, 16.4, 19.8,  0.0,  0.0, 10.0, 21.6, 15.0, 15.5, 
      14.0,  8.5, 11.5, 17.4, 15.6, 31.6, 20.6, 31.1, 16.3,  0.0, 33.1, 29.2,  11.2, 14.4, 60.0, 42.3,  9.5, 34.5, 
      36.5,  9.7,  0.0,  0.0,  7.6,  9.6,  9.3,  0.0)

#Northeastern
NE = c(25.3, 25.5, 24.0,  8.0,  0.0, 20.0,  0.0,  0.0,  0.0,  0.0, 33.8, 33.7, 15.1, 18.5, 44.8, 37.5,  
       0.0, 47.0, 20.0, 30.8, 30.0,  1.0,  1.2, 56.3, 0.0,  6.0,  0.0,  5.3,  7.2, 24.6, 30.0, 20.0, 14.5,  
       3.0, 28.0, 27.0, 50.0, 24.0, 30.0, 22.0, 16.0,  0.0,  0.0,  0.0,  0.0, 39.7,  9.3,  0.0,  2.1,  0.0, 
       46.7, 10.5,  0.0, 41.0, 10.3,  1.2, 23.9, 22.2, 24.1, 38.0,  9.0,  9.2,  6.6, 16.9, 10.0, 48.2,  6.5,  
       4.8, 25.0,  0.0,  0.0,  3.2,  44.0,  0.0,  9.3,  0.0, 20.0,  0.0,  4.8,  0.0,  0.0,  0.0,  0.0, 56.5, 
       39.2,  0.0,  6.4,  5.3,  0.0,  9.8,  0.0,  9.7,  4.5,  8.4,  0.8, 20.2,  0.0,  0.5,  5.3, 16.7, 45.2,  
       0.0,  0.6,  0.0,  3.1, 33.3,  6.0,  0.0, 13.2,  0.0, 21.0,  0.0,  8.4,  0.0,  0.0,  0.5,  4.5, 16.2,  
       0.0,  3.5,  20.0,  0.0,  0.0,  1.2,  0.0,  2.9,  0.0, 14.3,  0.0,  6.0,  0.0, 28.0,  0.0,  0.0,  0.0, 
       27.6, 33.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0, 15.0,   0.0,  0.0,  0.0, 25.7, 41.4, 41.6, 53.8, 48.5, 
       78.5, 12.7, 80.9, 67.0, 65.4, 21.0,  6.4,  0.0, 52.0, 45.0, 41.4, 14.3, 45.0, 27.0,  0.2, 30.0,  0.0,  0.0,  
       0.0,  0.0,  0.0,  1.0,  0.0,  0.0,  2.5,  0.0,  9.5,  0.0,  0.0,  0.0,  0.0,  2.1,  0.0,  0.0,  0.0,  0.0,  
       0.0,  7.2,  3.5,  0.0,  40.5, 25.8, 30.4,  0.0,  0.0,  0.0,  0.0, 0.0,  8.3,  0.0,  0.0,  0.0, 36.1,  0.0, 
       12.5,  0.0,  0.0,  0.0)

#Central 
C = c(39.6, 25.0,  0.0,  0.0, 29.7,  0.0,  3.1,  8.2,  3.2,  7.1,  0.0,  3.2,  4.2,  5.7, 30.0,  4.0,  0.0,  0.0,  
      0.0,  0.0,  0.0,  14.7,  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,  
      0.3,  0.5, 31.5,  8.2,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  0.0,  1.2,  0.0,  0.1,  1.0,  
      3.0,  0.5,  1.6)

#Eastern
E = c(0.0,   0.0,  26.5,  36.4,   0.0,   4.5,   0.0,   0.5,   0.7,  12.3,   0.5,   1.9,  66.4,  93.6,  68.7,  40.0,  65.0,  
      63.7,   9.2,   0.0,  0.0,  11.0,  69.6,  89.6, 160.0,  34.3,   0.0,  25.0,  19.5)

#Southern
S = c(27.9,   0.0,   3.4,   0.0,   0.8,  37.9,  32.4,  33.8,  15.8,   0.0,   0.0,   0.0,   0.0,  11.5,   1.7,   1.2,  21.2,   
      0.0,  30.0,   5.1,  2.5,   2.4,   5.0,   1.7,   0.0,   0.0,   2.1,   0.0,  10.5,  15.3,   0.0,   4.1,   9.0,  27.3,   
      6.5,   3.5,   0.0,   0.0,   0.0,   0.0,  3.6,   3.0,   1.0,   3.7,  15.6,  11.2,  24.0,   0.0,   0.0,  10.2,   0.0,   
      0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   7.0,  10.6,   3.6,   0.4,   3.8,   0.6,   0.0,  10.8,   5.0,  
      12.2,   3.6,   0.0,   8.8,   0.0,   0.0,   6.2,   0.0,  3.8,   0.0,   0.0,   0.0,  0.2,   0.0,   0.0,   0.4,   0.0,   
      0.0,   2.2,   0.0,   0.6,  76.6, 121.6,  60.0,   0.0, 114.6,   0.0,   0.0,   0.0,   0.0,  18.2,  40.4,  0.0,   0.0,  
      10.8,   0.0,   0.0,  35.0,   0.0,  33.5,  57.0,  10.5,   0.0,   0.0,   0.0,  30.8,  10.7,   0.0,  15.9,   0.0,   0.0)


##############################################################################
# R.code to construct histogram plots                                        #
##############################################################################
#Figure 5
png("Figure6.png", width = 6, height = 6, units = 'in', res = 300)
layout(matrix(c(1,1,2,2,3,3,0,4,4,5,5,0), 2, 6, byrow = TRUE))
hist(N, main="A", xlab = "Daily rainfall amounts",ylab = "Records",border = TRUE)
hist(NE, main="B", xlab = "Daily rainfall amounts",ylab = "Records")
hist(C, main="C", xlab = "Daily rainfall amounts",ylab = "Records")
hist(E, main="D", xlab = "Daily rainfall amounts",ylab = "Records")
hist(S, main="E", xlab = "Daily rainfall amounts",ylab = "Records")
dev.off() 

##############################################################################
# R.code to construct Normal Q-Q plots                                        #
##############################################################################
#Figure 7
png("Figure7.png", width = 6, height = 6, units = 'in', res = 300)
layout(matrix(c(1,1,2,2,3,3,0,4,4,5,5,0), 2, 6, byrow = TRUE))
log_N = log(N[N !=0])
log_NE = log(NE[NE !=0])
log_C = log(C[C !=0])
log_E = log(E[E !=0])
log_S = log(S[S !=0])

qqnorm(log_N, main = "A")
qqline(log_N,distribution = qnorm,probs = c(0.25, 0.75), qtype = 7)

qqnorm(log_NE, main = "B")
qqline(log_NE,distribution = qnorm, probs = c(0.25, 0.75), qtype = 7)

qqnorm(log_C, main = "C")
qqline(log_C,distribution = qnorm, probs = c(0.1, 0.85), qtype = 7)

qqnorm(log_E, main = "D")
qqline(log_E,distribution = qnorm, probs = c(0.20, 0.75), qtype = 7)

qqnorm(log_S, main = "E")
qqline(log_S,distribution = qnorm, probs = c(0.25, 0.75), qtype = 7)
dev.off()

##############################################################################
# R.code to estimate AIC results                                             #
##############################################################################

#Dataset: north
North_nonzero = N[N!=0]
cau_N = fitdistr(North_nonzero, "cauchy")
lognormal_N = fitdistr(North_nonzero, "lognormal")
normal_N = fitdistr(North_nonzero, "normal")
logistic_N = fitdistr(North_nonzero, "logistic")
t_N = fitdistr(North_nonzero, "t",df=length(North_nonzero)-1)

AIC(cau_N)
AIC(logistic_N)
AIC(lognormal_N)
AIC(normal_N)
AIC(t_N)

#Dataset: northeast
NE_nonzero = NE[NE !=0]
cau_NE = fitdistr(NE_nonzero, "cauchy")
lognormal_NE = fitdistr(NE_nonzero, "lognormal")
normal_NE = fitdistr(NE_nonzero, "normal")
logistic_NE = fitdistr(NE_nonzero, "logistic")
t_NE = fitdistr(NE_nonzero, "t",df=length(NE_nonzero)-1)

AIC(cau_NE)
AIC(logistic_NE)
AIC(lognormal_NE)
AIC(normal_NE)
AIC(t_NE)

#Dataset: central
C_nonzero = C[C !=0]
cau_C = fitdistr(C_nonzero, "cauchy")
lognormal_C = fitdistr(C_nonzero, "lognormal")
normal_C = fitdistr(C_nonzero, "normal")
logistic_C = fitdistr(C_nonzero, "logistic")
t_C = fitdistr(C_nonzero, "t",df=length(C_nonzero)-1)

AIC(cau_C)
AIC(logistic_C)
AIC(lognormal_C)
AIC(normal_C)
AIC(t_C)

#Dataset: east
E_nonzero = E[E!=0]
cau_E = fitdistr(E_nonzero, "cauchy")
lognormal_E = fitdistr(E_nonzero, "lognormal")
normal_E = fitdistr(E_nonzero, "normal")
logistic_E = fitdistr(E_nonzero, "logistic")
t_E = fitdistr(E_nonzero, "t", df=length(E_nonzero)-1)

AIC(cau_E)
AIC(logistic_E)
AIC(lognormal_E)
AIC(normal_E)
AIC(t_E)

#Dataset: south
S_nonzero = S[S!=0]
cau_S = fitdistr(S_nonzero, "cauchy")
lognormal_S = fitdistr(S_nonzero, "lognormal")
normal_S = fitdistr(S_nonzero, "normal")
logistic_S = fitdistr(S_nonzero, "logistic")
t_S = fitdistr(S_nonzero, "t",df=length(S_nonzero)-1)

AIC(cau_S)
AIC(logistic_S)
AIC(lognormal_S)
AIC(normal_S)
AIC(t_S)

AIC.N = data.frame()
AIC.NE = data.frame()
AIC.C = data.frame()
AIC.E = data.frame()
AIC.S = data.frame()
AIC.N = c(AIC(cau_N), AIC(logistic_N), AIC(lognormal_N), AIC(normal_N), AIC(t_N))
AIC.NE = c(AIC(cau_NE), AIC(logistic_NE), AIC(lognormal_NE), AIC(normal_NE), AIC(t_NE))
AIC.C = c(AIC(cau_C), AIC(logistic_C), AIC(lognormal_C), AIC(normal_C), AIC(t_C))
AIC.E = c(AIC(cau_E), AIC(logistic_E), AIC(lognormal_E), AIC(normal_E), AIC(t_E))
AIC.S = c(AIC(cau_S), AIC(logistic_S), AIC(lognormal_S), AIC(normal_S), AIC(t_S))
AIC = rbind(AIC.N, AIC.NE, AIC.C, AIC.E, AIC.S)

colnames(AIC) = c("Cauchy", "Logistic", "Lognormal", "Normal", "t")
AIC