### R-code for constructing FGCI, HPD Bayesian, and Standard Bootstrap ###

library("EnvStats")
library("HDInterval")
library("invgamma")
library("fitdistrplus")
library("MASS")
library("survival")
library("ggplot2")
library("magick")

CI.diffCVs <- function(M=15000,m=5000,k=3000,n1=25,n2=25,var1=0.5,var2=0.5,mean1=0,mean2=0,
                       delta1=0.5,delta2=0.5,alpha=0.05){
  starttime <- proc.time()
  
  cv1 <- sqrt(exp(var1)-1)                                  #variance of the LN distribution in population 1
  cv2 <- sqrt(exp(var2)-1)                                  #variance of the LN distribution in population 2
  ex1 <- delta1*exp(mean1+var1/2)                           #mean of the delta-LN distribution in population 1
  ex2 <- delta2*exp(mean2+var2/2)                           #mean of the delta-LN distribution in population 2
  var_x1 <- delta1*exp((2*mean1)+var1)*(exp(var1)-delta1)   #variance of the delta-LN distribution in population 1
  var_x2 <- delta2*exp((2*mean2)+var2)*(exp(var2)-delta2)   #variance of the delta-LN distribution in population 2
  eta1 <- sqrt(var_x1)/ex1                                  #cv of the delta-LN distribution in population 1
  eta2 <- sqrt(var_x2)/ex2                                  #cv of the delta-LN distribution in population 2
  gamma <- eta1-eta2                                        #difference between cvs of the delta-LN distribution
  
  #FGCI
  U1 <- matrix(rep(0,m))
  U2 <- matrix(rep(0,m))
  beta11 <- matrix(rep(0,m))
  beta12 <- matrix(rep(0,m))
  beta21 <- matrix(rep(0,m))
  beta22 <- matrix(rep(0,m))
  var1_GFQ <- matrix(rep(0,m))
  var2_GFQ <- matrix(rep(0,m))
  delta1_GFQ <- matrix(rep(0,m))
  delta2_GFQ <- matrix(rep(0,m))
  R1.GFQ <- matrix(rep(0,m))
  R2.GFQ <- matrix(rep(0,m))
  R.GFQ <- matrix(rep(0,m))
  
  #Bayesian 
  var1.linvj <- matrix(rep(0,m))
  var2.linvj <- matrix(rep(0,m))
  delta1.linvj <- matrix(rep(0,m))
  delta2.linvj <- matrix(rep(0,m))
  beta1.linvj <- matrix(rep(0,m))
  beta2.linvj <- matrix(rep(0,m))
  eta1.linvj <- matrix(rep(0,m))
  eta2.linvj <- matrix(rep(0,m))
  gamma.linvj <- matrix(rep(0,m))
  
  var1.jrule <- matrix(rep(0,m))
  var2.jrule <- matrix(rep(0,m))
  delta1.jrule <- matrix(rep(0,m))
  delta2.jrule <- matrix(rep(0,m))
  beta1.jrule <- matrix(rep(0,m))
  beta2.jrule <- matrix(rep(0,m))
  eta1.jrule <- matrix(rep(0,m))
  eta2.jrule <- matrix(rep(0,m))
  gamma.jrule <- matrix(rep(0,m))
  
  var1.uni <- matrix(rep(0,m))
  var2.uni <- matrix(rep(0,m))
  delta1.uni <- matrix(rep(0,m))
  delta2.uni <- matrix(rep(0,m))
  beta1.uni <- matrix(rep(0,m))
  beta2.uni <- matrix(rep(0,m))
  eta1.uni <- matrix(rep(0,m))
  eta2.uni <- matrix(rep(0,m))
  gamma.uni <- matrix(rep(0,m))
  
  #FGCI
  L.rFGCI <- numeric()
  U.rFGCI <- numeric ()
  CIr1 <- c(0,0)
  CP.rFGCI <- numeric()
  Length.rFGCI <- numeric()
  
  #HPD Bayesian using the left invariant Jeffreys prior
  hpd.linvj <- matrix(rep(0,2))
  L.rBaye_hpd.linvj <- numeric()
  U.rBaye_hpd.linvj <- numeric ()
  CIr2 <- c(0,0)
  CP.rBaye_hpd.linvj <- numeric()
  Length.rBaye_hpd.linvj <- numeric()
  
  #HPD Bayesian using the Jeffreys rule prior
  hpd.jrule <- matrix(rep(0,2))
  L.rBaye_hpd.jrule <- numeric()
  U.rBaye_hpd.jrule <- numeric ()
  CIr3 <- c(0,0)
  CP.rBaye_hpd.jrule <- numeric()
  Length.rBaye_hpd.jrule <- numeric()
  
  #HPD Bayesian using the uniform prior
  hpd.uni <- matrix(rep(0,2))
  L.rBaye_hpd.uni <- numeric()
  U.rBaye_hpd.uni <- numeric ()
  CIr4 <- c(0,0)
  CP.rBaye_hpd.uni <- numeric()
  Length.rBaye_hpd.uni <- numeric()
  
  #Standard Bootstrap
  L.rSB <- numeric()
  U.rSB <- numeric ()
  CIr5 <- c(0,0)
  CP.rSB <- numeric()
  Length.rSB <- numeric()
  
i = 0
  while(i < M){                                                #Begin loop i
    
    x1 <- rzmlnormAlt(n1,exp(mean1+var1/2),cv1,1-delta1)  
    x2 <- rzmlnormAlt(n2,exp(mean2+var2/2),cv2,1-delta2)
    
    nonzero_value1 <- log(x1[which(x1!=0)])                    #Find values of non-zero in population 1
    nonzero_value2 <- log(x2[which(x2!=0)])                    #Find values of non-zero in population 2
    
    size_nonzero1 <- length(nonzero_value1)                    #Find size of non-zero in population 1
    if(size_nonzero1 < n1*delta1){
      next 
    }
    
    size_nonzero2 <- length(nonzero_value2)                    #Find size of non-zero in population 2
    if(size_nonzero2 < n2*delta2){
      next 
    }
    
    mean_hat1 <- mean(nonzero_value1)
    var_hat1 <- var(nonzero_value1)
    delta_hat1 <- size_nonzero1/n1
    eta_hat1 <- sqrt((exp(var_hat1)-delta_hat1)/delta_hat1)
    
    mean_hat2 <- mean(nonzero_value2)
    var_hat2 <- var(nonzero_value2)
    delta_hat2 <- size_nonzero2/n2
    eta_hat2 <- sqrt((exp(var_hat2)-delta_hat2)/delta_hat2)
    
    #Fiducial generalized confidence interval (FGCI)...........................................................
    
    U1 <- rchisq(m,(size_nonzero1-1))                                   #Generate U1 from chi-square with n1(1)-1 df in population 1
    U2 <- rchisq(m,(size_nonzero2-1))                                   #Generate U2 from chi-square with n2(1)-1 df in population 2
    beta11 <- rbeta(m,size_nonzero1,(n1-size_nonzero1)+1)               #Generate beta11 from beta distribution in population 1  
    beta12 <- rbeta(m,size_nonzero1+1,n1-size_nonzero1)                 #Generate beta12 from beta distribution in population 2
    beta21 <- rbeta(m,size_nonzero2,(n2-size_nonzero2)+1)               #Generate beta21 from beta distribution in population 1
    beta22 <- rbeta(m,size_nonzero2+1,n2-size_nonzero2)                 #Generate beta22 from beta distribution in population 2
    var1_GFQ <- ((size_nonzero1-1)*var_hat1)/U1                         #Compute GFQ for delta-lognormal variance in population 1
    var2_GFQ <- ((size_nonzero2-1)*var_hat2)/U2                         #Compute GFQ for delta-lognormal variance in population 2
    delta1_GFQ <- (1/2)*beta11+(1/2)*beta12                             #Compute GFQ for probability of non-zero in population 1
    delta2_GFQ <- (1/2)*beta21+(1/2)*beta22                             #Compute GFQ for probability of non-zero in population 2
    R1.GFQ <- sqrt((exp(var1_GFQ)-delta1_GFQ)/delta1_GFQ)               #Compute GFQ for delta-lognormal coefficient of variation in population 1
    R2.GFQ <- sqrt((exp(var2_GFQ)-delta2_GFQ)/delta2_GFQ)               #Compute GFQ for delta-lognormal coefficient of variation in population 2
    R.GFQ <- R1.GFQ-R2.GFQ                                              #Compute GFQ for the difference cvs of delta-lognormal 
    if(is.na(R.GFQ)){
      next
    }
    L.rFGCI[i] <-quantile(R.GFQ,alpha/2)                                #Compute lower bounded based on FGCI
    U.rFGCI[i] <-quantile(R.GFQ,(1-alpha/2))                            #Compute upper bounded based on FGCI
    CIr1 <- rbind(c(L.rFGCI[i],U.rFGCI[i]))
    CP.rFGCI[i] <- ifelse(L.rFGCI[i]< gamma && gamma <U.rFGCI[i],1,0)   #Compute coverage probability
    Length.rFGCI[i] <- U.rFGCI[i]-L.rFGCI[i]                            #Compute lenght of interval
    
    #HPD Bayesian using the left invariant Jeffreys prior.......................................................

    beta1.linvj <- rbeta(m,size_nonzero1+(1/2),(n1-size_nonzero1)+(1/2))       #Generate the posterior distribution of having zero from beta distribution in population 1
    beta2.linvj <- rbeta(m,size_nonzero2+(1/2),(n2-size_nonzero2)+(1/2))       #Generate the posterior distribution of having zero from beta distribution in population 2
    var1.linvj <- rinvgamma(m, (size_nonzero1-1)/2,
                           ((size_nonzero1-1)*var_hat1)/2, scale = 1/rate)     #Generate the posterior distribution of variance from inverse gamma distribution in population 1
    var2.linvj <- rinvgamma(m, (size_nonzero2-1)/2,
                           ((size_nonzero2-1)*var_hat2)/2, scale = 1/rate)     #Generate the posterior distribution of variance from inverse gamma distribution in population 2
    delta1.linvj <- beta1.linvj                                                
    delta2.linvj <- beta2.linvj                                               
    eta1.linvj <-  sqrt((exp(var1.linvj)-delta1.linvj)/delta1.linvj)           #Compute delta-lognormal cv for group 1
    eta2.linvj <-  sqrt((exp(var2.linvj)-delta2.linvj)/delta2.linvj)           #Compute delta-lognormal cv for group 2
    gamma.linvj <- eta1.linvj-eta2.linvj
    hpd.linvj <- hdi(gamma.linvj,0.95)                                         #Compute hpd interval
    L.rBaye_hpd.linvj[i] <- hpd.linvj[1]                                       #Compute hpd lower bounded
    U.rBaye_hpd.linvj[i] <- hpd.linvj[2]                                       #Compute hpd upper bounded
    CIr2 <- rbind(c(L.rBaye_hpd.linvj[i],U.rBaye_hpd.linvj[i]))
    CP.rBaye_hpd.linvj[i] <- ifelse(L.rBaye_hpd.linvj[i]< gamma && gamma <U.rBaye_hpd.linvj[i],1,0)    #Compute coverage probability
    Length.rBaye_hpd.linvj[i] <- U.rBaye_hpd.linvj[i]-L.rBaye_hpd.linvj[i]                             #Compute lenght of interval
    
    
    #HPD Bayesian using the Jeffreys rule prior................................................................

    beta1.jrule <- rbeta(m,size_nonzero1+(3/2),(n1-size_nonzero1)+(1/2))        #Generate the posterior distribution of having zero from beta distribution in population 1
    beta2.jrule <- rbeta(m,size_nonzero2+(3/2),(n2-size_nonzero2)+(1/2))        #Generate the posterior distribution of having zero from beta distribution in population 2
    var1.jrule <- rinvgamma(m, size_nonzero1/2,
                            (size_nonzero1*var_hat1)/2, scale = 1/rate)         #Generate the posterior distribution of variance from inverse gamma distribution in population 1
    var2.jrule <- rinvgamma(m, size_nonzero2/2,
                            (size_nonzero2*var_hat2)/2, scale = 1/rate)         #Generate the posterior distribution of variance from inverse gamma distribution in population 2
    delta1.jrule <- beta1.jrule                                                 
    delta2.jrule <- beta2.jrule                                                 
    eta1.jrule <-  sqrt((exp(var1.jrule)-delta1.jrule)/delta1.jrule)            #Compute delta-lognormal cv for group 1
    eta2.jrule <-  sqrt((exp(var2.jrule)-delta2.jrule)/delta2.jrule)            #Compute delta-lognormal cv for group 2
    gamma.jrule <- eta1.jrule-eta2.jrule
    hpd.jrule <- hdi(gamma.jrule,0.95)                                          #Compute hpd interval
    L.rBaye_hpd.jrule[i] <- hpd.jrule[1]                                        #Compute hpd lower bounded
    U.rBaye_hpd.jrule[i] <- hpd.jrule[2]                                        #Compute hpd upper bounded 
    CIr3 <- rbind(c(L.rBaye_hpd.jrule[i],U.rBaye_hpd.jrule[i]))
    CP.rBaye_hpd.jrule[i] <- ifelse(L.rBaye_hpd.jrule[i]< gamma && gamma <U.rBaye_hpd.jrule[i],1,0)     #Compute coverage probability
    Length.rBaye_hpd.jrule[i] <- U.rBaye_hpd.jrule[i]-L.rBaye_hpd.jrule[i]                              #Compute lenght of interval  
    
    #HPD Bayesian using the uniform prior.......................................

    beta1.uni <- rbeta(m,size_nonzero1+1,(n1-size_nonzero1)+1)                  #Generate the posterior distribution of having zero from beta distribution in population 1
    beta2.uni <- rbeta(m,size_nonzero2+1,(n2-size_nonzero2)+1)                  #Generate the posterior distribution of having zero from beta distribution in population 2
    var1.uni <- rinvgamma(m, (size_nonzero1-2)/2,
                          ((size_nonzero1-2)*var_hat1)/2, scale = 1/rate)       #Generate the posterior distribution of variance from inverse gamma distribution in population 1
    var2.uni <- rinvgamma(m, (size_nonzero2-2)/2,
                          ((size_nonzero2-2)*var_hat2)/2, scale = 1/rate)       #Generate the posterior distribution of variance from inverse gamma distribution in population 2
    delta1.uni <- beta1.uni                                                           
    delta2.uni <- beta2.uni                                                           
    eta1.uni <-  sqrt((exp(var1.uni)-delta1.uni)/delta1.uni)                    #Compute delta-lognormal cv for group 1
    eta2.uni <-  sqrt((exp(var2.uni)-delta2.uni)/delta2.uni)                    #Compute delta-lognormal cv for group 2
    gamma.uni <- eta1.uni-eta2.uni
    hpd.uni <- hdi(gamma.uni,0.95)                                              #Compute hpd interval
    L.rBaye_hpd.uni[i] <- hpd.uni[1]                                            #Compute lower bounded
    U.rBaye_hpd.uni[i] <- hpd.uni[2]                                            #Compute upper bounded
    CIr4 <- rbind(c(L.rBaye_hpd.uni[i],U.rBaye_hpd.uni[i]))
    CP.rBaye_hpd.uni[i] <- ifelse(L.rBaye_hpd.uni[i]< gamma && gamma <U.rBaye_hpd.uni[i],1,0)       #Compute coverage probability
    Length.rBaye_hpd.uni[i] <- U.rBaye_hpd.uni[i]-L.rBaye_hpd.uni[i]                                #Compute lenght of interval
    
    #Standard bootstrap.........................................................................................
    
    y1y2 <- data.frame(cbind(nonzero_value1,nonzero_value2))
    y1.B <- with(y1y2,matrix(sample(nonzero_value1[1:size_nonzero1],size = size_nonzero1*k,replace = TRUE),k,size_nonzero1))   #Resampling sample for group 1
    y1.B.var <- apply(y1.B,1,var)
    y1.B.cv <- sqrt((exp(y1.B.var)-delta_hat1)/delta_hat1)          #Compute cv in population 1
    y2.B <- with(y1y2,matrix(sample(nonzero_value2,size = size_nonzero2*k,replace = TRUE),k,size_nonzero2))                    #Resampling sample for group 2
    y2.B.var <- apply(y2.B,1,var)
    y2.B.cv <- sqrt((exp(y2.B.var)-delta_hat2)/delta_hat2)          #Compute cv in population 2
    B.diff <- y1.B.cv-y2.B.cv                                       #Compute the difference cvs
    CI.SB <- (eta_hat1-eta_hat2)+c(-1,1)*qnorm(0.975,0,1)*sd(B.diff)
    L.rSB[i] <- CI.SB[1]                                            #Compute lower bounded
    U.rSB[i] <- CI.SB[2]                                            #Compute upper bounded
    CIr5 <- rbind(c(L.rSB[i],U.rSB[i]))
    CP.rSB[i] <- ifelse(L.rSB[i]< gamma && gamma <U.rSB[i],1,0)     #Compute coverage probability
    Length.rSB[i] <- U.rSB[i]-L.rSB[i]                              #Compute length of interval
    
    i=i+1
    print(i)
    
  }                                                                     #end loop i
  
  ACP.rFGCI <- mean(CP.rFGCI)	                                          #compute coverage probability for FGCI
  Ex.length.rFGCI <- mean(Length.rFGCI)                                 #compute average length for FGCI
  ACP.rBaye_hpd.linvj <- mean(CP.rBaye_hpd.linvj)	                      #compute coverage probability for HPD Bayesian (left invariant Jeffreys prior)
  Ex.length.rBaye_hpd.linvj <- mean(Length.rBaye_hpd.linvj)             #compute average length for HPD Bayesian (left invariant Jeffreys prior)
  ACP.rBaye_hpd.jrule <- mean(CP.rBaye_hpd.jrule)	                      #compute coverage probability for HPD Bayesian (Jeffreys' Rule prior)
  Ex.length.rBaye_hpd.jrule <- mean(Length.rBaye_hpd.jrule)             #compute average length for HPD Bayesian (Jeffreys' Rule prior)
  ACP.rBaye_hpd.uni <- mean(CP.rBaye_hpd.uni)	                          #compute coverage probability for HPD Bayesian (uniform prior)
  Ex.length.rBaye_hpd.uni <- mean(Length.rBaye_hpd.uni)                 #compute average length for HPD Bayesian (uniform prior)
  ACP.rSB <- mean(CP.rSB)	                                              #compute coverage probability based on SB
  Ex.Length.rSB <- mean(Length.rSB)                                   	#compute average length based on SB
  elapseTime = proc.time()-starttime
  cat("M =",M,"m =",m,"k =",k,"n1 =",n1,"n2 =",n2,"delta1 =",delta1,"delta2 =",delta2,"var1 =",var1,"var2 =",var2,"\n",
      "E1(x_0)=",n1*delta1,"E2(x_0)=",n2*delta2,"cv1=",cv1,"cv2=",cv2,"gamma",gamma,"\n",
      "CP|rFGCI =",ACP.rFGCI,"Length|rFGCI =",Ex.length.rFGCI,"\n",
      "CP|rBaye_hpd.linvj =",ACP.rBaye_hpd.linvj,"Length|rBaye_hpd.linvj =",Ex.length.rBaye_hpd.linvj,"\n",
      "CP|rBaye_hpd.jrule =",ACP.rBaye_hpd.jrule,"Length|rBaye_hpd.jrule =",Ex.length.rBaye_hpd.jrule,"\n",
      "CP|rBaye_hpd.uni =",ACP.rBaye_hpd.uni,"Length|rBaye_hpd.uni =",Ex.length.rBaye_hpd.uni,"\n",
      "CP|rSB =",ACP.rSB,"Length|rSB =",Ex.Length.rSB)
  result <- c("M =",M,"m =",m,"k =",k,"n1 =",n1,"n2 =",n2,"delta1 =",delta1,"delta2 =",delta2,"var1 =",var1,"var2 =",var2,"\n",
              "E1(x_0)=",n1*delta1,"E2(x_0)=",n2*delta2,"cv1=",cv1,"cv2=",cv2,"gamma=",gamma,"\n",
              "CP|rFGCI =",ACP.rFGCI,"Length|rFGCI =",Ex.length.rFGCI,"\n",
              "CP|rBaye_hpd.linvj =",ACP.rBaye_hpd.linvj,"Length|rBaye_hpd.linvj =",Ex.length.rBaye_hpd.linvj,"\n",
              "CP|rBaye_hpd.jrule =",ACP.rBaye_hpd.jrule,"Length|rBaye_hpd.jrule =",Ex.length.rBaye_hpd.jrule,"\n",
              "CP|rBaye_hpd.uni =",ACP.rBaye_hpd.uni,"Length|rBaye_hpd.uni =",Ex.length.rBaye_hpd.uni,"\n",
              "CP|rSB =",ACP.rSB,"Length|rSB =",Ex.Length.rSB,"\n",
              "Time =",elapseTime["elapsed"])
  return(result)  
}

### R-code for sets of real data (rainfall series) ###

# Histogram....................................................................................................

Jana <- c(0.5,11.5,15.2,0,3.5,85.3,47.1,62,27,49.9,713.9,267.8,25.6,0,89,43.2,96.6,15.5,85.1,9.5,18,30.7,678.7,41.4,
          24.6,0,12.5,0,6,19,38.7,0,84,297.3,388,372.7,156.7,0,136.3,40.5,15,87.9,106.5,151.7,79.4,101.6,444.1,640.5,
          290.2,10,49.2,181.1,32,54,7,78.5,147.6,154.6,104.7,272,62.9,152,0,14.7,129.5,80.1,95.4,103.1,50.5,162.7,581.3,
          324.1,12,0,0,5.5,52.9,81.4,25,199.2,60,443.8,429.5,573,8,0,0,20,55,27,29.1,95,145.5,185,247,49.2,92.5,28,0,0,
          39,30,24.6,70.9,70.8,127.1,211,656.7,352.4,5.5,51.9,254.5,136.8,52.9,62,205,153.3,165.1,1017.3,269.9)
Ranod <- c(177.2,37,0,66.7,60.5,28.1,15.6,79.8,5.4,89,856.6,223.3,32.2,5,92,81,47.5,0,8.5,6.8,0,42,295,25.1,0,0,0,0,0,0,
           8,0,8,314.8,470.5,381.2,85.3,0,193.5,26.7,0,0,6.7,6,16.7,311.1,518.9,406.9,428.3,0,37.5,5.5,5.5,26.4,0,32,39.9,
           216.5,216.8,279.1,56.8,152.5,0,107.4,157.5,77.3,18.8,62.7,8.7,251.9,630.9,291.1,1.6,0,0,0,62.2,0,0,36,18.4,197.5,
           478.8,473.5,0,0,0,38.2,41.5,0,25,12.5,76.1,279.7,652.8,125,0,40.7,0,0,21.5,16.2,7.2,5,0,44,186.4,1324.5,830.6,
           26.3,20.7,29.4,23.1,0,0,1.2,105,54,1059.1,345.6)

Jana_nonzero <- Jana[which(Jana!=0)]
Ranod_nonzero <- Ranod[which(Ranod!=0)]

#Histogram, theoretical densities, and AIC values for positive monthly rainfall data in Jana and Ranod.........

#Jana
fn_Jana <- fitdist(Jana_nonzero, "norm")                        #Normal distribution
fca_Jana <- fitdist(Jana_nonzero, "cauchy")                     #Cauchy distribution
fln_Jana <- fitdist(Jana_nonzero, "lnorm")                      #Lognormal distribution
fexp_Jana <- fitdist(Jana_nonzero, "exp")                       #Exponential distribution
fgam_Jana <- fitdist(Jana_nonzero, "gamma")                     #Gamma distribution
summary(fn_Jana)
summary(fca_Jana)
summary(fln_Jana)
summary(fexp_Jana)
summary(fgam_Jana)
plot.legend_Jana <- c("normal", "cauchy", "lognormal", "exponential", "gamma")
png("Fig4a.png", width = 1200, height = 1200, units = "px", pointsize = 12, res = 300)
denscomp(list(fn_Jana, fca_Jana, fln_Jana, fexp_Jana, fgam_Jana), xlim=c(0,1200), 
         ylim=c(0,0.012), probability = TRUE, main="Densities of positive rainfall data \n from Jana",
         xlab="Rainfall (mm.)", legendtext = plot.legend_Jana)
mtext('A', side=1, line = 4, cex = 1.2, at = 9)
dev.off()

#Ranod
fn_Ranod <- fitdist(Ranod_nonzero, "norm")
fca_Ranod <- fitdist(Ranod_nonzero, "cauchy")
fln_Ranod <- fitdist(Ranod_nonzero, "lnorm")
fexp_Ranod <- fitdist(Ranod_nonzero, "exp")
fgam_Ranod <- fitdist(Ranod_nonzero, "gamma")
summary(fn_Ranod)
summary(fca_Ranod)
summary(fln_Ranod)
summary(fexp_Ranod)
summary(fgam_Ranod)
plot.legend_Ranod <- c("normal", "cauchy", "lognormal", "exponential", "gamma")
png("Fig4b.png", width = 1200, height = 1200, units = "px", pointsize = 12, res = 300)
denscomp(list(fn_Ranod, fca_Ranod, fln_Ranod, fexp_Ranod, fgam_Ranod), xlim=c(0,1400), 
         ylim=c(0,0.012), probability = TRUE, main = "Densities of positive rainfall data \n from Ranod",
         xlab="Rainfall (mm.)", legendtext = plot.legend_Ranod)
mtext('B', side=1, line = 4, cex = 1.2, at = 9)
dev.off()

#Normal Q-Q plots of log-transformed of positive monthly rainfall data........................................

#Jana
png("Fig5a.png", width = 1200, height = 1200, units = "px", pointsize = 12, res = 300)
qqnorm(log(Jana_nonzero), main = "Normal Q-Q plot of log-transformed \n positive rainfall data from Jana", xlab = "Theoretical Quantiles", ylab = "Sample Quantiles", plot.it = TRUE,
       datax = FALSE)
qqline(log(Jana_nonzero), datax = FALSE, distribution = qnorm, probs = c(0.25, 0.75), qtype = 7)
mtext('A', side=1, line = 4, cex = 1.2, at = -2.5)
dev.off()

#Ranod
png("Fig5b.png", width = 1200, height = 1200, units = "px", pointsize = 12, res = 300)
qqnorm(log(Ranod_nonzero), main = "Normal Q-Q plot of log-transformed \n positive rainfall data from Ranod", xlab = "Theoretical Quantiles", ylab = "Sample Quantiles", plot.it = TRUE,
       datax = FALSE)
qqline(log(Ranod_nonzero), datax = FALSE, distribution = qnorm, probs = c(0.25, 0.75), qtype = 7)
mtext('B', side=1, line = 4, cex = 1.2, at = -2.5)
dev.off()

#Histogram and theoretical densities for log-transformed of positive monthly rainfall data....................

#Jana
fn_Jana.log <- fitdist(log(Jana_nonzero), "norm")
png("Fig6a.png", width = 1200, height = 1200, units = "px", pointsize = 12, res = 300)
denscomp(fn_Jana.log, xlim=c(0,8), ylim=c(0,0.4), probability = TRUE, main="Density of log-transformed \n positive rainfall data from Jana",
         xlab="Rainfall (mm.)", addlegend = FALSE)
mtext('A', side=1, line = 4, cex = 1.2, at = 0)
dev.off()

#Ranod
fn_Ranod.log <- fitdist(log(Ranod_nonzero), "norm")
png("Fig6b.png", width = 1200, height = 1200, units = "px", pointsize = 12, res = 300)
denscomp(fn_Ranod.log, xlim=c(0,8), ylim=c(0,0.3), probability = TRUE, main="Density of log-transformed \n positive rainfall data from Ranod",
         xlab="Rainfall (mm.)", addlegend = FALSE)
mtext('B', side=1, line = 4, cex = 1.2, at = 0)
dev.off()

# confidence intervals and Bayesian credible intervals plotted...............................................

Data <- data.frame(Methods=c("FGCI","B_linvj","B_jrule","B_uni","SB"),
                   diff.cvs=c(-0.0716, -0.0716, -0.0716, -0.0716, -0.0716),
                   lower=c(-4.0492, -3.7924, -3.5602, -3.7008, -2.9163),
                   upper=c(-0.0558, 0.1359, 0.1619, 0.0646, -0.1118))
head(Data)

png("Fig7.png", width = 1200, height = 1200, units = "px", pointsize = 12, res = 300)
ggplot(data = Data, aes(x = Methods, y = diff.cvs)) +
  geom_point() +
  geom_errorbar(aes(x = Methods, ymin = lower, ymax = upper),
                color="black",width=0.5,size=0.5) +
  ggtitle("Plot of confidence intervals \n and Bayesian credible intervals") +
  theme_bw() +
  theme(plot.title = element_text(hjust = 0.5)) +
  labs(x= "Methods", y = "Confidence/Bayesian credible intervals")
dev.off()

### Combine images ###

F4A = image_read("Fig4a.png")
F4B = image_read("Fig4b.png")
Fig4 = c(F4A, F4B)
side_by_side = image_append(Fig4, stack=F)
image_write(side_by_side, path = "Fig4.png", format = "png")

F5A = image_read("Fig5a.png")
F5B = image_read("Fig5b.png")
Fig5 = c(F5A, F5B)
side_by_side = image_append(Fig5, stack=F)
image_write(side_by_side, path = "Fig5.png", format = "png")

F6A = image_read("Fig6a.png")
F6B = image_read("Fig6b.png")
Fig6 = c(F6A, F6B)
side_by_side = image_append(Fig6, stack=F)
image_write(side_by_side, path = "Fig6.png", format = "png")