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

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

CI.ratioVarDLN = function(M=5000, m=2500, n1=50, n2=50, mu1=3, mu2=3, sigmaSq1=1, 
                          sigmaSq2=1, delta1=0.1, delta2=0.1, alpha=0.05){    
  
  # 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_nonzero1 = At least the sample size of nonzero in Group1 = n1-(n1*delta1)
  # Temp.size_nonzero2 = At least the sample size of nonzero in Group2 = n2-(n2*delta2)
  # var_deltaLN        = variance of DLN
  # ratio_VarDLN       = logarithm of ratio of delta-lognormal variances
  
  # mu_GPQ            = generalized pivotal quantity (GPQ) of mu
  # sigmaSq_GPQ       = GPQ of sigmaSq
  # delta_GPQ         = GPQ of probability of zero values
  # R.VST             = GPQ of log-ratio of delta-lognormal variances
  # L.GCI             = lower limit of generalized confidence interval (GCI)
  # U.GCI             = upper limit of GCI
  # CP.GCI            = coverage probability based on GCI
  # Length.GCI        = length of GCI
  # ACP.GCI           = average coverage probability of GCI	      		                       
  # Ex.length.GCI     = average length of GCI
  # RAL.GCI           = relative average length of GCI
  
  # mu_FGPQ           = fiducial generalised pivotal quantity (FGPQ) of mu
  # sigmaSq_FGPQ      = FGPQ of sigmaSq
  # deltanon_FGPQ     = FGPQ of probability of non-zero values
  # T_var             = FGPQ of variance of DLN
  # T_ratio.FGCI      = FGPQ of log-ratio of delta-lognormal variances
  # 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
  # RAL.FGCI          = relative average length of FGCI
  # 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
  
  # HPD.Jef          = highest posterior density interval based on Jeffreys' prior (HPD-Jef)
  # post_sigmaSq     = posterior of sigmaSq based on Jeffreys' prior
  # post_mu          = posterior of mu based on Jeffreys' prior
  # post_delta       = posterior of delta based on Jeffreys' prior
  # post_varDLN      = posterior of delta-lognormal variance
  # ratio.post.Jef   = ratio of two posteriors of delta-lognormal variance based on Jeffreys' prior
  # L.HPD.Jef        = lower limit of HPD-Jef
  # U.HPD.Jef        = upper limit of HPD-Jef
  # CP.HPD.Jef       = coverage probability of HPD-Jef
  # Length.HPD.Jef   = average length of HPD-Jef
  # ACP.HPD.Jef      = average coverage probability of HPD-Jef      		                       
  # Ex.length.HPD.Jef = average length of HPD-Jef
  
  # HPD.rul             = highest posterior density interval based on Jeffreys' Rule prior (HPD-Rul)
  # post_sigmaSq.rul    = posterior of sigmaSq based on Jeffreys' Rule prior
  # post_mu.rul         = posterior of mu based on Jeffreys' Rule prior
  # post_delta.rul      = posterior of delta based on Jeffreys' Rule prior
  # post_varDLN.rul     = posterior of delta-lognormal variance
  # ratio.post.rul      = ratio of two posteriors of delta-lognormal variance based on Jeffreys' Rule prior
  # L.HPD.rul           = lower limit of HPD-Rul
  # U.HPD.rul           = upper limit of HPD-Rul
  # CP.HPD.rul          = coverage probability of HPD-Rul
  # Length.HPD.rul      = average length of HPD-Rul
  # ACP.HPD.rul         = average coverage probability of HPD-Rul      		                       
  # Ex.length.HPD.rul   = average length of HPD-Rul
  
  # HPD.NG              = highest posterior density interval based on normal-gamma prior (HPD-NG)
  # post_sigmaSq.NG     = posterior of sigmaSq based on Jeffreys' Rule prior
  # post_mu.NG          = posterior of mu based on Jeffreys' Rule prior
  # post_deltanon.NG    = posterior of delta based on Jeffreys' Rule prior
  # post_varDLN.NG      = posterior of delta-lognormal variance
  # ratio.post.NG       = ratio of two posteriors of delta-lognormal variance
  # L.HPD.NG            = lower limit of HPD-NG
  # U.HPD.NG            = upper limit of HPD-NG
  # CP.HPD.NG           = coverage probability of HPD-NG
  # Length.HPD.NG       = average length of HPD-NG
  # ACP.HPD.NG          = average coverage probability of HPD-NG      		                       
  # Ex.length.HPD.NG    = average length of HPD-NG
  
  Ex_LN1         = exp(mu1+(sigmaSq1/2))
  cv_LN1          = sqrt(exp(sigmaSq1)-1)                                                
  var_deltaLN1 = (1-delta1)*exp((2*mu1)+sigmaSq1)*(exp(sigmaSq1)+delta1-1)
  
  Ex_LN2         = exp(mu2+(sigmaSq2/2))
  cv_LN2          = sqrt(exp(sigmaSq2)-1)                                                 
  var_deltaLN2 = (1-delta2)*exp((2*mu2)+sigmaSq2)*(exp(sigmaSq2)+delta2-1)
  
  ratio_VarDLN = log(var_deltaLN1/var_deltaLN2)
  
  # GCI
  U1.VST = matrix(rep(0,m))
  U2.VST = matrix(rep(0,m))
  Z1.VST = matrix(rep(0,m))
  Z2.VST = matrix(rep(0,m))
  V1.VST = matrix(rep(0,m))
  V2.VST = matrix(rep(0,m))
  mu_GPQ1 = matrix(rep(0,m))
  mu_GPQ2 = matrix(rep(0,m))
  sigmaSq_GPQ1 = matrix(rep(0,m))
  sigmaSq_GPQ2 = matrix(rep(0,m))
  delta_GPQ1 = matrix(rep(0,m))
  delta_GPQ2 = matrix(rep(0,m))
  R.VST = matrix(rep(0,m))
  L.GCI = numeric()
  U.GCI = numeric()
  CP.GCI = numeric()
  Length.GCI = numeric()
  
  # FGCI
  mu_FGPQ1 = matrix(rep(0,m))
  mu_FGPQ2 = matrix(rep(0,m))
  sigmaSq_FGPQ1 = matrix(rep(0,m))
  sigmaSq_FGPQ2 = matrix(rep(0,m))
  delta_FGPQ1 = matrix(rep(0,m))
  delta_FGPQ2 = matrix(rep(0,m))
  T_var1 = matrix(rep(0,m))
  T_var2 = matrix(rep(0,m))
  T_ratio.FGCI = matrix(rep(0,m))
  L.FGCI = numeric()
  U.FGCI = numeric()
  CP.FGCI = numeric()
  Length.FGCI = numeric()
  
  # MOVER
  L.mover = numeric()
  U.mover = numeric()
  CP.mover = numeric()
  Length.mover = numeric()
  
  # HPD-Jef
  post_mu1 = rep(0,m)
  post_sigmaSq1= rep(0,m)
  post_delta1 = rep(0,m)
  post_varDLN1 = rep(0,m)
  post_mu2 = rep(0,m)
  post_sigmaSq2= rep(0,m)
  post_delta2 = rep(0,m)
  post_varDLN2 = rep(0,m)
  ratio.post.Jef = rep(0,m)
  L.HPD.Jef = numeric()
  U.HPD.Jef = numeric()
  CP.HPD.Jef = numeric()
  Length.HPD.Jef = numeric()
  
  # HPD-Rul
  post_mu.rul1 = rep(0,m)
  post_sigmaSq.rul1 = rep(0,m)
  post_delta.rul1  = rep(0,m)
  post_varDLN.rul1  = rep(0,m)
  
  post_mu.rul2 = rep(0,m)
  post_sigmaSq.rul2 = rep(0,m)
  post_delta.rul2  = rep(0,m)
  post_varDLN.rul2  = rep(0,m)
  L.HPD.rul  = numeric()
  U.HPD.rul  = numeric()
  CP.HPD.rul  = numeric()
  Length.HPD.rul  = numeric()
  
  
  # HPD-NG
  post_mu.NG1 = rep(0,m)
  post_sigmaSq.NG1 = rep(0,m)
  post_deltanon.NG1  = rep(0,m)
  post_mu.NG2 = rep(0,m)
  post_sigmaSq.NG2 = rep(0,m)
  post_deltanon.NG2  = rep(0,m)
  L.HPD.NG  = numeric()
  U.HPD.NG  = numeric()
  CP.HPD.NG  = numeric()
  Length.HPD.NG  = numeric()
  
  i = 1
  while(i <= M){                                                   #begin loop i
    
    x1 = rzmlnormAlt(n1, Ex_LN1, cv_LN1, delta1)                   #generate x1 from DLN in Population 1
    x2 = rzmlnormAlt(n2, Ex_LN2, cv_LN1, delta2)                   #generate x2 from DLN in Population 2
    
    Temp.nonzero_value1 = log(x1[which(x1!=0)])                    #find values is not zero in Population 1
    Temp.nonzero_value2 = log(x2[which(x2!=0)])                    #find values is not zero in Population 2
    
    
    Temp.size_nonzero1 = length(Temp.nonzero_value1)  
    if(Temp.size_nonzero1 < (n1-(n1*delta1))){next}
    
    Temp.size_nonzero2 = length(Temp.nonzero_value2)  
    if(Temp.size_nonzero2 < (n2-(n2*delta2))){next}
    
    nonzero_value1 = Temp.nonzero_value1
    size_nonzero1 = Temp.size_nonzero1
    size_zero1 = n1-size_nonzero1
    sum_nonzero1 = sum(nonzero_value1^2)
    mu_hat1 = mean(nonzero_value1)
    sigmaSq_hat1 =var(nonzero_value1)
    delta_hat1 = 1-(size_nonzero1/n1)
    varDLN_hat1 = (1-delta_hat1)*exp(2*mu_hat1+sigmaSq_hat1)*(exp(sigmaSq_hat1)-(1-delta_hat1))
    
    nonzero_value2 = Temp.nonzero_value2
    size_nonzero2 = Temp.size_nonzero2
    size_zero2 = n2-size_nonzero2
    sum_nonzero2 = sum(nonzero_value2^2)
    mu_hat2 = mean(nonzero_value2)
    sigmaSq_hat2 =var(nonzero_value2)
    delta_hat2 = 1-(size_nonzero2/n2)
    varDLN_hat2 = (1-delta_hat2)*exp(2*mu_hat2+sigmaSq_hat2)*(exp(sigmaSq_hat2)-(1-delta_hat2))
    
    ##############################################################################
    #                         Algorithm 1: GCI                                   #
    ##############################################################################
    
    U1.VST = rchisq(m,(size_nonzero1-1))                     #generate U1 from chi-square with n1(1)-1 df in Population 1
    Z1.VST = rnorm(m,0,1)                                    #generate Z1 from standard normal in Population 1
    V1.VST = rnorm(m,0,1)                                    #generate w3 from standard normal
    
    U2.VST = rchisq(m,(size_nonzero2-1))                     #generate U2 from chi-square with n2(1)-1 df in Population 2
    Z2.VST = rnorm(m,0,1)                                    #generate Z2 from standard normal in Population 2
    V2.VST = rnorm(m,0,1)                                    #generate w3 from standard normal
    
    mu_GPQ1 = mu_hat1-(V1.VST*sqrt(((size_nonzero1-1)*sigmaSq_hat1)/(size_nonzero1*U1.VST)))                      
    sigmaSq_GPQ1 = ((size_nonzero1-1)*sigmaSq_hat1)/U1.VST                                                                
    delta_GPQ1 = (sin(asin(sqrt(delta_hat1))-(Z1.VST/(2*sqrt(n1)))))^2                                               
    R1.VST =(1-delta_GPQ1)*exp(2*mu_GPQ1+sigmaSq_GPQ1)*(exp(sigmaSq_GPQ1)-(1-delta_GPQ1))              
    
    mu_GPQ2 = mu_hat2-(V2.VST*sqrt(((size_nonzero2-1)*sigmaSq_hat2)/(size_nonzero2*U2.VST)))
    sigmaSq_GPQ2 = ((size_nonzero2-1)*sigmaSq_hat2)/U2.VST
    delta_GPQ2 = (sin(asin(sqrt(delta_hat2))-(Z2.VST/(2*sqrt(n2)))))^2                                                       
    R2.VST =(1-delta_GPQ2)*exp(2*mu_GPQ2+sigmaSq_GPQ2)*(exp(sigmaSq_GPQ2)-(1-delta_GPQ2))
    R.VST = log(R1.VST/R2.VST)                                               
    
    L.GCI[i] =quantile(R.VST,0.025)                                       
    U.GCI[i] =quantile(R.VST,0.975)                                      
    CP.GCI[i] = ifelse(L.GCI[i]< ratio_VarDLN && ratio_VarDLN <U.GCI[i],1,0)           
    Length.GCI[i] = U.GCI[i]-L.GCI[i]                                  
    
    ##############################################################################
    #                         Algorithm 2: FGCI                                  #
    ##############################################################################
    
    z1 = rnorm(m,0,1) 
    u1 = rchisq(m,size_nonzero1-1)/(size_nonzero1-1)
    UU1 = sqrt(u1)
    mu_FGPQ1 =  mu_hat1-(z1*sqrt(sigmaSq_hat1)/(UU1*sqrt(size_nonzero1)))
    sigmaSq_FGPQ1 =  sigmaSq_hat1/u1
    delta_FGPQ1 = rbeta(m,size_nonzero1+0.5, size_zero1+0.5)
    T_var1 = delta_FGPQ1*exp((2*mu_FGPQ1 )+sigmaSq_FGPQ1 )*(exp(sigmaSq_FGPQ1)-delta_FGPQ1) 
    
    z2 = rnorm(m,0,1) 
    u2 = rchisq(m,size_nonzero2-1)/(size_nonzero2-1)
    UU2 = sqrt(u2)
    mu_FGPQ2 =  mu_hat2-(z2*sqrt(sigmaSq_hat2)/(UU2*sqrt(size_nonzero2)))
    sigmaSq_FGPQ2 =  sigmaSq_hat2/u2
    delta_FGPQ2 = rbeta(m,size_nonzero2+0.5, size_zero2+0.5)
    T_var2 = delta_FGPQ2*exp((2*mu_FGPQ2 )+sigmaSq_FGPQ2 )*(exp(sigmaSq_FGPQ2)-delta_FGPQ2) 
    T_ratio.FGCI = log(T_var1/T_var2) 
    
    L.FGCI[i] =quantile(T_ratio.FGCI,0.025)                                       
    U.FGCI[i] =quantile(T_ratio.FGCI,0.975)                                       
    CP.FGCI[i] = ifelse(L.FGCI[i]< ratio_VarDLN && ratio_VarDLN <U.FGCI[i],1,0)
    Length.FGCI[i] = U.FGCI[i]-L.FGCI[i]   
    
    ##############################################################################
    #                          Algorithm 3:  MOVER                               #
    ##############################################################################
    
    #Group1
    #CIs for ln(1-delta)
    z1 = qnorm(1-(alpha/2),0,1)
    delta_hatnon1 = 1-delta_hat1
    delta1.1 = delta_hatnon1+(0.5*z1^2/n1)
    delta2.1 = sqrt((delta_hatnon1*(1-delta_hatnon1)+ (0.25*z1^2/n1))/n1)
    delta3.1 = 1+(z1^2/n1) 
    g_hat1.1 = log(delta_hatnon1)
    g.1l.1 = log((delta1.1-(z1*delta2.1))/delta3.1)
    g.1u.1 = log((delta1.1+(z1*delta2.1))/delta3.1)
    
    #CIs for 2mu+sigmaSq
    chi.l.1 = qchisq(1-(alpha/2),size_nonzero1-1)
    chi.u.1 = qchisq(alpha/2,size_nonzero1-1)
    g_hat2.1 = (2*mu_hat1)+sigmaSq_hat1
    mean.Nor.1 = sigmaSq_hat1*z1^2/size_nonzero1
    sigmaSq.Nor.l.1 =  (sigmaSq_hat1^2)*(1-((size_nonzero1-1)/chi.l.1))^2
    sigmaSq.Nor.u.1 =  (sigmaSq_hat1^2)*(((size_nonzero1-1)/chi.u.1)-1)^2
    g.2l.1 = g_hat2.1 - sqrt(4*mean.Nor.1+sigmaSq.Nor.l.1)
    g.2u.1 = g_hat2.1 + sqrt(4*mean.Nor.1+sigmaSq.Nor.u.1)
    
    #CIs for ln(exp(sigmaSq)-(1-delta))
    g_hat3.1 = exp(sigmaSq_hat1)- delta_hatnon1
    sd_delta_hatnon1 = z1^2*((delta_hatnon1*(1-delta_hatnon1)+ (0.25*z1^2/n1))/n1)/(1+(z1^2/n1))^2
    g.3l.1 = log(g_hat3.1 - sqrt(exp(sigmaSq.Nor.l.1)+sd_delta_hatnon1))
    g.3u.1 = log(g_hat3.1 + sqrt(exp(sigmaSq.Nor.u.1)+sd_delta_hatnon1))
    
    #Group1: CIs for ln(1-delta)+(2mu+sigmaSq)+ln(exp(sigmaSq)-(1-delta))
    l1 = (g_hat1.1+g_hat2.1+log(g_hat3.1)) - sqrt((g_hat1.1-g.1l.1)^2+(g_hat2.1-g.2l.1)^2+(log(g_hat3.1)-g.3l.1)^2) 
    u1 = (g_hat1.1+g_hat2.1+log(g_hat3.1)) + sqrt((g_hat1.1-g.1u.1)^2+(g_hat2.1-g.2u.1)^2+(log(g_hat3.1)-g.3u.1)^2) 
    
    #Group2
    #CIs for ln(1-delta)
    z2 = qnorm(1-(alpha/2),0,1)
    delta_hatnon2 = 1-delta_hat2
    delta1.2 = delta_hatnon2+(0.5*z2^2/n2)
    delta2.2 = sqrt((delta_hatnon2*(1-delta_hatnon2)+ (0.25*z2^2/n2))/n2)
    delta3.2 = 1+(z2^2/n2) 
    g_hat1.2 = log(delta_hatnon2)
    g.1l.2 = log((delta1.2-(z2*delta2.2))/delta3.2)
    g.1u.2 = log((delta1.2+(z2*delta2.2))/delta3.2)
    
    #CIs for 2mu+sigmaSq
    chi.l.2 = qchisq(1-(alpha/2),size_nonzero2-1)
    chi.u.2 = qchisq(alpha/2,size_nonzero2-1)
    g_hat2.2 = (2*mu_hat2)+sigmaSq_hat2
    mean.Nor.2 = sigmaSq_hat2*z2^2/size_nonzero2
    sigmaSq.Nor.l.2 =  (sigmaSq_hat2^2)*(1-((size_nonzero2-1)/chi.l.2))^2
    sigmaSq.Nor.u.2 =  (sigmaSq_hat2^2)*(((size_nonzero2-1)/chi.u.2)-1)^2
    g.2l.2 = g_hat2.2 - sqrt(4*mean.Nor.2+sigmaSq.Nor.l.2)
    g.2u.2 = g_hat2.2 + sqrt(4*mean.Nor.2+sigmaSq.Nor.u.2)
    
    #CIs for ln(exp(sigmaSq)-(1-delta))
    g_hat3.2 = exp(sigmaSq_hat2)-delta_hatnon2
    sd_delta_hatnon2 = z2^2*((delta_hatnon2*(1-delta_hatnon2)+ (0.25*z2^2/n2))/n2)/(1+(z2^2/n2))^2
    g.3l.2 = log(g_hat3.2 - sqrt(exp(sigmaSq.Nor.l.2)+sd_delta_hatnon2))
    g.3u.2 = log(g_hat3.2 + sqrt(exp(sigmaSq.Nor.u.2)+sd_delta_hatnon2))
    
    #Group2: CIs for ln(1-delta)+(2mu+sigmaSq)+ln(exp(sigmaSq)-(1-delta))
    l2 = (g_hat1.2+g_hat2.2+log(g_hat3.2)) - sqrt((g_hat1.2-g.1l.2)^2+(g_hat2.2-g.2l.2)^2+(log(g_hat3.2)-g.3l.2)^2) 
    u2 = (g_hat1.2+g_hat2.2+log(g_hat3.2)) + sqrt((g_hat1.2-g.1u.2)^2+(g_hat2.2-g.2u.2)^2+(log(g_hat3.2)-g.3u.2)^2) 
    
    v_hat1 = log(1-delta_hat1)+(2*mu_hat1+sigmaSq_hat1)+log(exp(sigmaSq_hat1)-(1-delta_hat1))
    v_hat2 = log(1-delta_hat2)+(2*mu_hat2+sigmaSq_hat2)+log(exp(sigmaSq_hat2)-(1-delta_hat2))
    
    L.mover[i] =  (v_hat1-v_hat2) - sqrt((v_hat1-l1)^2+(u2-v_hat2)^2)
    U.mover[i] =  (v_hat1-v_hat2) + sqrt((u1-v_hat1)^2+(v_hat2-l2)^2)
    CP.mover[i] = ifelse(L.mover[i] < ratio_VarDLN && ratio_VarDLN < U.mover[i],1,0)
    Length.mover[i] = U.mover[i] - L.mover[i]
    if(is.na(Length.mover[i])|is.infinite(Length.mover[i])){next}
    
    
    ##############################################################################
    #      Algorithm 4: HPD intervals based on different priors                 #
    ##############################################################################
    
    ############################### HPD-Jef ######################################    
    
    v1 = size_nonzero1-1
    shap1 = v1/2
    scal1 = v1*sigmaSq_hat1/2
    post_sigmaSq1 = rinvgamma(m, shap1, scal1)
    Tem.post_mu1 = numeric()
    for(k in 1:length(post_sigmaSq1)){
      Tem.post_mu1[k] = rnorm(1,mu_hat1,sqrt(post_sigmaSq1[k]/size_nonzero1))
    }
    post_mu1 = Tem.post_mu1
    post_delta1 = rbeta(m, size_zero1+0.5, size_nonzero1+0.5)   
    post_varDLN1 = (1-post_delta1)*exp(2*post_mu1+post_sigmaSq1)*(exp(post_sigmaSq1)-(1-post_delta1))
    
    v2 = size_nonzero2-1
    shap2 = v2/2
    scal2 = v2*sigmaSq_hat2/2
    post_sigmaSq2 = rinvgamma(m, shap2, scal2)
    Tem.post_mu2 = numeric()
    for(k in 1:length(post_sigmaSq2)){
      Tem.post_mu2[k] = rnorm(1,mu_hat2,sqrt(post_sigmaSq2[k]/size_nonzero2))
    }
    post_mu2 = Tem.post_mu2
    post_delta2 = rbeta(m, size_zero2+0.5, size_nonzero2+0.5)   
    post_varDLN2 = (1-post_delta2)*exp(2*post_mu2+post_sigmaSq2)*(exp(post_sigmaSq2)-(1-post_delta2))
    
    ratio.post.Jef = log(post_varDLN1/post_varDLN2)
    
    HPD.Jef = hdi(ratio.post.Jef, credMass = 1-alpha , allowSplit=FALSE)
    L.HPD.Jef[i] = HPD.Jef[1]
    U.HPD.Jef[i] = HPD.Jef[2]
    CP.HPD.Jef[i] = ifelse(L.HPD.Jef[i]< ratio_VarDLN && ratio_VarDLN < U.HPD.Jef[i],1,0)
    Length.HPD.Jef[i] = U.HPD.Jef[i]-L.HPD.Jef[i]
    
    ############################### HPD-Rul ######################################    
    
    v1 = size_nonzero1+1
    shap1 = v1/2
    scal1 = v1*sigmaSq_hat1/2
    post_sigmaSq.rul1 = rinvgamma(m, shap1, scal1)
    Tem.post_mu.rul1 = numeric()
    for(k in 1:length(post_sigmaSq.rul1)){
      Tem.post_mu.rul1[k] = rnorm(1,mu_hat1,sqrt(post_sigmaSq.rul1[k]/size_nonzero1))
    }
    post_mu.rul1 = Tem.post_mu.rul1
    post_delta.rul1  = rbeta(m, size_zero1+0.5, size_nonzero1+1.5)   
    post_varDLN.rul1  = (1-post_delta.rul1)*exp(2*post_mu.rul1+post_sigmaSq.rul1)*(exp(post_sigmaSq.rul1)-(1-post_delta.rul1))
    
    v2 = size_nonzero2+1
    shap2 = v2/2
    scal2 = v2*sigmaSq_hat2/2
    post_sigmaSq.rul2 = rinvgamma(m, shap2, scal2)
    Tem.post_mu.rul2 = numeric()
    for(k in 1:length(post_sigmaSq.rul2)){
      Tem.post_mu.rul2[k] = rnorm(1,mu_hat2,sqrt(post_sigmaSq.rul2[k]/size_nonzero2))
    }
    post_mu.rul2 = Tem.post_mu.rul2
    post_delta.rul2  = rbeta(m, size_zero2+0.5, size_nonzero2+1.5)   
    post_varDLN.rul2  = (1-post_delta.rul2)*exp(2*post_mu.rul2+post_sigmaSq.rul2)*(exp(post_sigmaSq.rul2)-(1-post_delta.rul2))
    
    ratio.post.rul = log(post_varDLN.rul1/post_varDLN.rul2)
    if(sum(is.na(ratio.post.rul)*1)){next}
    
    HPD.rul = hdi(ratio.post.rul, credMass = 1-alpha ,allowSplit=FALSE)
    L.HPD.rul[i] = HPD.rul[1]
    U.HPD.rul[i] = HPD.rul[2]
    CP.HPD.rul[i] = ifelse(L.HPD.rul[i]< ratio_VarDLN && ratio_VarDLN < U.HPD.rul[i],1,0)
    Length.HPD.rul[i] = U.HPD.rul[i]-L.HPD.rul[i]
    
    ############################### HPD-NG ######################################    
    
    v1 = size_nonzero1-1
    alpha_n1 = (v1-1)/2
    beta_n1 = (sum_nonzero1-(size_nonzero1*mu_hat1^2))/2
    post_sigmaSq.NG1 = rinvgamma(m, alpha_n1, beta_n1)
    post_mu.NG1 = rst(1,mu_hat1,sigmaSq_hat1/size_nonzero1,v1)
    k_alpha2 = (qnorm(0.025,0,1))^2
    post_deltanon.NG1  = rbeta(m,size_nonzero1+(1/3)+(k_alpha2/6),size_zero1+(1/3)+(k_alpha2/6))
    post_varDLN.NG1  = post_deltanon.NG1*exp(2*post_mu.NG1+post_sigmaSq.NG1)*(exp(post_sigmaSq.NG1)-post_deltanon.NG1)
    
    v2 = size_nonzero2-1
    alpha_n2 = (v2-1)/2
    beta_n2 = (sum_nonzero2-(size_nonzero2*mu_hat2^2))/2
    post_sigmaSq.NG2 = rinvgamma(m, alpha_n2, beta_n2)
    post_mu.NG2 = rst(1,mu_hat2,sigmaSq_hat2/size_nonzero2,v2)
    post_deltanon.NG2  = rbeta(m,size_nonzero2+(1/3)+(k_alpha2/6), size_zero2+(1/3)+(k_alpha2/6))
    post_varDLN.NG2  = post_deltanon.NG2*exp(2*post_mu.NG2+post_sigmaSq.NG2)*(exp(post_sigmaSq.NG2)-post_deltanon.NG2)
    ratio.post.NG = log(post_varDLN.NG1/post_varDLN.NG2)
    
    HPD.NG  = hdi(ratio.post.NG, credMass = 1-alpha , allowSplit=FALSE)
    L.HPD.NG[i] = HPD.NG[1]
    U.HPD.NG[i] = HPD.NG[2]
    CP.HPD.NG[i] = ifelse(L.HPD.NG[i]< ratio_VarDLN && ratio_VarDLN < U.HPD.NG[i],1,0)
    Length.HPD.NG[i] = U.HPD.NG[i]-L.HPD.NG[i]
    
    # Count after overall runs is achieved
    i=i+1
    if(i%%1000==0){
      print(i)
    }
  }                                                                          						 #end loop i
  
  
  ACP.HPD.rul = mean(CP.HPD.rul)*100	      		                       
  Ex.length.HPD.rul = mean(Length.HPD.rul) 
  
  ACP.HPD.Jef = mean(CP.HPD.Jef)*100	      		                       
  Ex.length.HPD.Jef = mean(Length.HPD.Jef) 		                       
  RAL.length.HPD.Jef = mean(Length.HPD.Jef)/Ex.length.HPD.rul  	
  
  ACP.HPD.NG= mean(CP.HPD.NG)*100	      		                           
  Ex.length.HPD.NG = mean(Length.HPD.NG)   		                 
  RAL.length.HPD.NG = Ex.length.HPD.NG/Ex.length.HPD.rul  
  
  ACP.GCI = mean(CP.GCI)*100	                                         
  Ex.length.GCI = mean(Length.GCI)                                    
  RAL.length.GCI = Ex.length.GCI/Ex.length.HPD.rul                       
  
  ACP.FGCI = mean(CP.FGCI)*100                                         
  Ex.length.FGCI = mean(Length.FGCI)                                   
  RAL.length.FGCI = Ex.length.FGCI/Ex.length.HPD.rul                   
  
  ACP.mover = mean(CP.mover)*100	                                     
  Ex.length.mover = mean(Length.mover)                                
  RAL.length.mover = Ex.length.mover/Ex.length.HPD.rul                  
  
  cat("M=",M, "m =",m, "n=",c(n1, n2),"\n",
      "mean",c(mu1, mu2),"\n",
      "sigmaSq=",c(sigmaSq1, sigmaSq2),"\n",
      "delta=",c(delta1, delta2),"\n",
      "ratio_VarDLN=", ratio_VarDLN,"\n", 
      "CP|HPD.Jef=", ACP.HPD.Jef, "RAL.HPD.Jef=", RAL.length.HPD.Jef,"\n",
      "CP|HPD.rul=", ACP.HPD.rul,"\n",
      "CP|HPD.NG=", ACP.HPD.NG, "RAL.HPD.NG=", RAL.length.HPD.NG,"\n",
      "CP|GCI =", ACP.GCI, "RAL.GCI =", RAL.length.GCI,"\n",
      "CP|FGCI =", ACP.FGCI, "RAL.FGCI =", RAL.length.FGCI,"\n",
      "CP|mover =", ACP.mover, "RAL.mover =", RAL.length.mover)
}

##############################################################################
# R.code to construct the normal Q-Q plots of rainfall data                  #
##############################################################################
# north region
N  = c(5.0, 11.7,  0.0,  0.0,  0.0, 88.8, 43.0, 17.2,  7.6,  0.0,  0.0,  6.2,  0.0,  0.0,  9.7,  7.4,  
       7.4,  0.0,  0.0,  1.2,  0.0, 0.0, 14.3,  7.1,  5.3, 16.2, 14.5,  2.4, 30.0,  0.0,  2.8,  0.0,  
       0.0, 10.0, 0.0, 4.8, 21.0, 0.0, 30.0, 18.2, 2.8, 25.9, 77.0, 0.0, 27.2, 0.0, 30.1, 17.6,  
       0.0, 12.1, 42.4, 0.0, 14.5, 25.3, 22.1, 3.4, 0.0, 9.4, 20.1, 20.2, 19.8, 0.0)

# northeast region
NE = c(50.8,  20.0,  19.4,  14.0,   0.0,   0.0,  15.2,  40.0,  12.5,   0.0,   4.5,   0.0,   3.3,   6.7,   
       0.4,   6.4,   3.4,  7.8,   2.5,   0.7,  10.6,   0.0,   4.4,   5.4,   5.4,   0.0,   4.7,  30.0,  
       27.0,   0.0,  31.0,   0.0,   7.0,  11.0,  10.0,   8.0,  25.0,   2.4,  30.0,   0.0,   7.0,  21.2,   
       0.0,   4.0,  30.7,  21.5,  26.8,   5.1,   9.7,  27.0,  33.2,  43.2,  12.8,   0.0,   8.0,   4.6,   
       3.5,   1.4,   0.5,  15.2,  99.5,   5.9, 128.8, 149.2,  78.7,  68.0,  67.7,  48.3,  156.5,  64.5,  
       98.0,   5.3,  10.8,   0.0,   0.0,   5.2,  12.3,   1.3,  11.0,   7.0,  13.0,   0.0,  20.0,  19.2,   
       3.6,  0.0,   1.1,   0.0,   9.3,   0.0,  13.5,   0.0,   4.2,  60.5,  68.1,  65.2,  29.0,  33.1,  
       65.4,  83.0, 135.0,   0.0,  7.4,   0.0,   8.2,   0.0,   0.0,   2.0,   0.0,   0.0,   8.0,   0.0,   
       0.0,   0.0,  14.5,  30.3,  31.5,  11.5,  20.0,  14.5,   4.0,  10.0,   8.0,   0.0,   0.0,   1.4,   
       6.4,   2.3,   0.0,   0.0,   1.5,  17.8,   0.0,  17.8,   0.0,  62.9,  17.4,   0.0,   0.0,   0.0,   
       0.0,   0.0,   0.0,   0.0,  38.0,   4.2,  10.0,  23.1,  76.1,  12.0,  25.5,  45.0,  28.7,  0.0,  
       49.5,  11.7,  38.3,   0.0,  53.4, 141.9,  23.8,  41.6,  76.1,  10.5,   0.0,   0.0,  57.5,   0.5,   
       2.8,  35.1,  0.0,  12.5,  22.0,   8.4,   4.5,  37.0,  17.5,  10.3,  18.8,   8.0,   0.0,   0.0,   
       0.0,   7.2,   0.0,  60.0,  14.2,  0.0,   1.0,   3.4,  25.5,   0.0,   4.5,  11.2,   6.5,   0.0,   
       0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,   0.0,  0.0,   6.4,   0.0,   0.0,   0.0,   0.0)

nonzero_value_N = log(N[which(N!=0)])
nonzero_value_NE = log(NE[which(NE!=0)])

png("Fig1.png", width = 8, height = 4, units = 'in', res = 300)
par(mfrow=c(1,2))
qqnorm(nonzero_value_N, main = "A", cex.main=1.5, cex.lab=1.5, cex.axis=1.5)
qqline(nonzero_value_N,distribution = qnorm, probs = c(0.15, 0.85), qtype = 6)

qqnorm(nonzero_value_NE, main = "B", cex.main=1.5, cex.lab=1.5, cex.axis=1.5)
qqline(nonzero_value_NE,distribution = qnorm, probs = c(0.05, 0.90), qtype = 6)
dev.off()

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

png("Fig2.png", width = 8, height = 4, units = 'in', res = 300)
par(mfrow=c(1,2))
hist(nonzero_value_N, main="A", xlab = "Log-positive rainfall amounts", ylab = "Records", 
     breaks = 5, cex.main=1.5, cex.lab=1.5, cex.axis=1.5)
hist(nonzero_value_NE, main="B", xlab = "Log-positive rainfall amounts", ylab = "Records", 
     breaks = 6, cex.main=1.5, cex.lab=1.5, cex.axis=1.5)
dev.off()  


##############################################################################
# R.code to estimate AIC results                                             #
##############################################################################
a1 <- N[which(N!=0)]
expo1 <- fitdistr(a1,"exponential")
weibull1 = fitdistr(a1,"weibull")
lognorm1 <- fitdistr(a1,"lognormal")
norm1 <- fitdistr(a1,"normal")
t1 <- fitdistr(a1,"t")
cauchy1 <- fitdistr(a1,"cauchy")
logis1 <- fitdistr(a1,"logistic")

a2 <- NE[which(NE!=0)]
expo2 <- fitdistr(a2,"exponential")
weibull2 = fitdistr(a2,"weibull")
lognorm2 <- fitdistr(a2,"lognormal")
norm2 <- fitdistr(a2,"normal")
t2 <- fitdistr(a2,"t")
cauchy2 <- fitdistr(a2,"cauchy")
logis2 <- fitdistr(a2,"logistic")

AIC.N = data.frame()
AIC.NE = data.frame()
AIC.N = c(AIC(expo1), AIC(weibull1), AIC(lognorm1), AIC(norm1),
          AIC(t1), AIC(cauchy1), AIC(logis1))
AIC.NE = c(AIC(expo2), AIC(weibull2),AIC(lognorm2), AIC(norm2),
           AIC(t2), AIC(cauchy2), AIC(logis2))
AIC = rbind(AIC.N,AIC.NE)
colnames(AIC) = c("Exponential", "Weibull", "Lognormal", "Normal", "t", "Cauchy", "Logistic")
AIC