##############################################################################
# R code to apply with rainfall data 
# For the Difference Mean of Delta-gamma Distributions (Unequal Sample Sizes)
##############################################################################

library(stats)
library(invgamma)
library(EnvStats)
library(HDInterval)
library(extraDistr)
library(MASS)
library(fitdistrplus)
library(LaplacesDemon)

# Difference


# x = data
# y = data
# M = number of simulation runs
# m = number of generalized computation  
# n_x0 = size of zero in population x
# n_y0 = size of zero in population y
# n_x1 = size of non-zero in population x
# n_y1 = size of non-zero in population y
# n_x = number of sample size of x 
# n_y = number of sample size of y 
# shape_x = shape parameter of x
# shape_y = shape parameter of y
# rate_x = rate parameter of x
# rate_y = rate parameter of y
# theta = difference between two means
# alpha = level of significance 

# L.CI.J = lower CI based on credible confidence interval Jeffreys prior
# U.CI.J = upper CI based on credible confidence interval Jeffreys prior
# L.CI.HPD.J = lower CI based on HPD Jeffreys prior
# U.CI.HPD.J = upper CI based on HPD Jeffreys prior
# L.CI.U = lower CI based on credible confidence interval uniform prior
# U.CI.U = upper CI based on credible confidence interval uniform prior
# L.CI.HPD.U = lower CI based on HPD uniform prior
# U.CI.HPD.U = upper CI based on HPD uniform prior
# L.CI.F = lower CI based on fiducial quantities
# U.CI.F = upper CI based on fiducial quantities

set.seed(1) 
m	<-  5000

beta.jef_x  <- rep(0,m)
ex.jef_x    <- rep(0,m)
var.jef_x   <- rep(0,m)
delta.jef_x <- rep(0,m)
beta.jef_y  <- rep(0,m)
ex.jef_y    <- rep(0,m)
var.jef_y   <- rep(0,m)
delta.jef_y <- rep(0,m)
theta.jef   	<- rep(0,m)

L.CI.J 		<- numeric()
U.CI.J 		<- numeric()
CP.J 		<- numeric()
Length.J 	<- numeric()

HPD.J	 	<- matrix(rep(0,2))
L.CI.HPD.J 	<- numeric()
U.CI.HPD.J 	<- numeric()
CP.HPD.J 	<- numeric()
Length.HPD.J <- numeric()

beta.uni_x  <- rep(0,m)
ex.uni_x    <- rep(0,m)
var.uni_x   <- rep(0,m)
delta.uni_x <- rep(0,m)
beta.uni_y  <- rep(0,m)
ex.uni_y    <- rep(0,m)
var.uni_y   <- rep(0,m)
delta.uni_y <- rep(0,m)
theta.uni   <- rep(0,m)

L.CI.U 		<- numeric()
U.CI.U 		<- numeric()
CP.U 		<- numeric()
Length.U 	<- numeric()

HPD.U	 	<- matrix(rep(0,2))
L.CI.HPD.U 	<- numeric()
U.CI.HPD.U 	<- numeric()
CP.HPD.U 	<- numeric()
Length.HPD.U <- numeric()

beta0_x 	<- rep(0,m)
beta1_x		<- rep(0,m)
R.delta_x	<- rep(0,m)
R.Fmu_x		<- rep(0,m)
R.Fsigsq_x	<- rep(0,m)
R.mean_x	<- rep(0,m)
beta0_y 	<- rep(0,m)
beta1_y		<- rep(0,m)
R.delta_y	<- rep(0,m)
R.Fmu_y		<- rep(0,m)
R.Fsigsq_y	<- rep(0,m)
R.mean_y	<- rep(0,m)
R.theta 	<- rep(0,m)

L.CI.F		<- numeric()
U.CI.F		<- numeric()
CP.F		<- numeric()
Length.F 	<- numeric()

x_MLE	 	<- matrix(rep(0,2))
y_MLE	 	<- matrix(rep(0,2))

##############################################################################
# Start Iteration 
##############################################################################

#Dataset 4

D4 =	c(0,0,0,22.5,0,0,3.5,6,13.2,0
		,51.1,29.9,34.3,0,42.5,0,2.5,0,0,0
		,46.8,0,0,0,0,28.7,0,0,7.3,0
		,0,0,0,16,4.7,11,0,5.2,24.5,61.2
		,61.5,7.4,34,0,102.2,0,0,57.2,0,0,0)

  x 	<- D4^(1/3)
  n_x 	<- length(D4)
  
  nonzero_value_x 	<-  x[which(x!=0)] 
  n_nonzero_x 		<- length(nonzero_value_x)
  n_zero_x 			<- n_x - n_nonzero_x
  
  x_MLE 			<- fitdistr(nonzero_value_x, "gamma", start=list(shape=1, rate=1))$estimate
  
  shape_x 			<- x_MLE[1]						
  rate_x			<- x_MLE[2]				
  scale_x   		<- 1/rate_x
  
#Dataset 5

D5 =	c(0.4,0.1,0.3,1.7,3.3,0,0,0,0,0
		,88.6,0,59.2,34.5,0,0,0,0,0,0
		,7.4,0,0,0,121.5,2.8,0,0,0.5,0
		,0,0,0,0,0,0,7.5,0,25.6,2
		,0,0,0,11,0,0,7.7,12.1,38,0
		,55.5,28,30.4,0,48.4,0,6,13,8.3,0
		,0,0,0,0,0.1,0,0,0,11.4,1.6
		,54.5,0,0,0,0,0,38.3,0,0,0
		,0,0,13.6,0,0,0,0,0,31.7,0
		,0,39,13.2,0,0,0,0.5,0,0,0
		,11,0,0,0,0,27.1,0,0,3.9,0
		,0,0,0,7.4,13.4,0,3.2,1.1,0,43
		,0,12.6,40.4,19.2,0,5.8,29.8,0,0,9.1
		,0,0,0,0.5,35,0,0,0,2.7,9.1
		,0,0,0.8,136.2,0.5,21.6,38.5,0,19.5,44.5
		,83.5,0,24.5,0,29,33.2,1.5,0,26.3,20.2
		,61.4,11.8,71,0,29.1,0,0,15.6,0,0,0.9)

  y 	<- D5^(1/3)
  n_y 	<- length(D5)
  
  nonzero_value_y 	<-  y[which(y!=0)] 
  n_nonzero_y 		<- length(nonzero_value_y)
  n_zero_y 			<- n_y - n_nonzero_y
  
  y_MLE 			<- fitdistr(nonzero_value_y, "gamma", start=list(shape=1, rate=1))$estimate
  
  shape_y 			<- y_MLE[1]						
  rate_y			<- y_MLE[2]				
  scale_y   		<- 1/rate_y
  
  delta_x 	<- n_zero_x/(n_zero_x+n_nonzero_x)
  delta_y 	<- n_zero_y/(n_zero_y+n_nonzero_y)

  theta_x		<- (1-delta_x)*shape_x/rate_x
  theta_y		<- (1-delta_y)*shape_y/rate_y
  
  theta  		<- theta_x - theta_y
  alpha			<- 0.05
  
  mu_x 		<- (shape_x/scale_x)^1/3*(1-(1/(9*shape_x)))
  sigsq_x	<-  1/(9*(shape_x)^(1/3)*(scale_x^(2/3)))
  
  mu_y		<- (shape_y/scale_y)^1/3*(1-(1/(9*shape_y)))
  sigsq_y	<-  1/(9*(shape_y)^(1/3)*(scale_y^(2/3)))
 
  xbar 		<- mean(nonzero_value_x)
  s_x 		<- sd(nonzero_value_x) 
  sum_x    	<- sum((nonzero_value_x-mean(nonzero_value_x))^2)
  ex_x_hat 	<- xbar
  var_x_hat <- (s_x)^2
  u_x 		<- rchisq(m , n_nonzero_x-1)/(n_nonzero_x-1)
  z_x 		<- rnorm(m)
  Fmu_x 	<- xbar + z_x*s_x/sqrt(n_nonzero_x)/sqrt(u_x)
  Fsigsq_x 	<- (s_x)^2/u_x
  
  ybar 		<- mean(nonzero_value_y)
  s_y 		<- sd(nonzero_value_y) 
  sum_y    	<- sum((nonzero_value_y-mean(nonzero_value_y))^2) 
  ex_y_hat 	<- ybar
  var_y_hat <- (s_y)^2
  u_y 		<- rchisq(m , n_nonzero_y-1)/(n_nonzero_y-1)
  z_y 		<- rnorm(m)
  Fmu_y 	<- ybar + z_y*s_y/sqrt(n_nonzero_y)/sqrt(u_y)
  Fsigsq_y 	<- (s_y)^2/u_y
  
##############################################################################
# compute CI based on credible confidence interval Jeffreys prior
##############################################################################

beta.jef_x  <- rbeta(m,(n_x -n_nonzero_x )+(1/2),n_nonzero_x +(3/2)) 
ex.jef_x    <- rnorm(m,ex_x_hat,sigsq_x/n_nonzero_x)
var.jef_x  	<- rinvgamma(m, n_nonzero_x/2,sum_x/2)	 
delta.jef_x <- beta.jef_x  
mean.jef_x  <- (0.5*ex.jef_x  + sqrt(0.25*(ex.jef_x )^2 + var.jef_x ))^3 
  
beta.jef_y  <- rbeta(m,(n_y -n_nonzero_y )+(1/2),n_nonzero_y +(3/2))   
ex.jef_y    <- rnorm(m,ex_y_hat,sigsq_y/n_nonzero_y)
var.jef_y  	<- rinvgamma(m, n_nonzero_y/2,sum_y/2)   
delta.jef_y <- beta.jef_y  
mean.jef_y  <- (0.5*ex.jef_y  + sqrt(0.25*(ex.jef_y )^2 + var.jef_y ))^3 

theta.jef   <- (delta.jef_x * mean.jef_x) - (delta.jef_y * mean.jef_y)

L.CI.J  	<- quantile(theta.jef,alpha/2)                             
U.CI.J  	<- quantile(theta.jef,(1-alpha/2)) 

Length.J	<- U.CI.J - L.CI.J

##############################################################################
# compute CI based on HPD Jeffreys prior
##############################################################################

HPD.J	 		<- hdi(theta.jef,0.95)                                                
L.CI.HPD.J 		<- HPD.J[1]                                            	     
U.CI.HPD.J 		<- HPD.J[2]                                            
        
Length.HPD.J	<- U.CI.HPD.J - L.CI.HPD.J
  
##############################################################################
# compute CI based on credible confidence interval uniform prior
##############################################################################
 
beta.uni_x 		<- rbeta(m,(n_x-n_nonzero_x)+1, n_nonzero_x+1) 
ex.uni_x   		<- rnorm(m,ex_x_hat,sigsq_x/n_nonzero_x)
var.uni_x 		<- rinvgamma(m, (n_nonzero_x-3)/2,sum_x/2)
delta.uni_x 	<- beta.uni_x                                                           		
mean.uni_x  	<- (0.5*ex.uni_x + sqrt(0.25*(ex.uni_x)^2 + var.uni_x))^3

beta.uni_y 		<- rbeta(m,(n_y-n_nonzero_y)+1, n_nonzero_y+1)  
ex.uni_y   		<- rnorm(m,ex_y_hat,sigsq_y/n_nonzero_y)	
var.uni_y 		<- rinvgamma(m, (n_nonzero_y-3)/2,sum_y/2)	
delta.uni_y 	<- beta.uni_y                                                           		
mean.uni_y  	<- (0.5*ex.uni_y + sqrt(0.25*(ex.uni_y)^2 + var.uni_y))^3  
   
theta.uni		<- (delta.uni_x * mean.uni_x) - (delta.uni_y * mean.uni_y)
  
L.CI.U  		<- quantile(theta.uni,alpha/2)                                  	
U.CI.U	 		<- quantile(theta.uni,(1-alpha/2))  
     
Length.U	 	<- U.CI.U - L.CI.U
 
##############################################################################
# compute CI based on HPD uniform prior
##############################################################################

HPD.U 			<- hdi(theta.uni,0.95)                                                    	
L.CI.HPD.U		<- HPD.U[1]                                             
U.CI.HPD.U		<- HPD.U[2]                                               

Length.HPD.U	<- U.CI.HPD.U - L.CI.HPD.U	
  
##############################################################################
# compute CI based on fiducial quantities
##############################################################################

beta0_x 		<- rbeta(m,n_nonzero_x ,n_zero_x +1)
beta1_x 		<- rbeta(m,n_nonzero_x +1,n_zero_x)
R.delta_x 		<- (1/2)*beta0_x  + (1/2)*beta1_x 
R.Fmu_x 		<- xbar + (z_x*s_x) /sqrt(n_nonzero_x )/sqrt(u_x )
R.Fsigsq_x 		<- (s_x)^2/u_x
R.mean_x 		<- (0.5*R.Fmu_x  + sqrt(0.25*(R.Fmu_x )^2 + R.Fsigsq_x ))^3
  
beta0_y 		<- rbeta(m,n_nonzero_y ,n_zero_y +1)
beta1_y 		<- rbeta(m,n_nonzero_y +1,n_zero_y)
R.delta_y 		<- (1/2)*beta0_y  + (1/2)*beta1_y 
R.Fmu_y 		<- ybar + (z_y*s_y) /sqrt(n_nonzero_y )/sqrt(u_y )
R.Fsigsq_y 		<- (s_y)^2/u_y
R.mean_y 		<- (0.5*R.Fmu_y  + sqrt(0.25*(R.Fmu_y )^2 + R.Fsigsq_y ))^3
  
R.theta 		<- (R.delta_x * R.mean_x) - (R.delta_y * R.mean_y)

L.CI.F  		<- quantile(R.theta,alpha/2,type=8,na.rm = TRUE)       			
U.CI.F 			<- quantile(R.theta,(1-alpha/2),type=8,na.rm = TRUE)			

Length.F	 	<- U.CI.F - L.CI.F
	
##############################################################################

c(n_x, n_y, delta_x, delta_y)

c(shape_x, rate_x)

c(shape_y, rate_y)

c(theta, theta_x, theta_y)

c(L.CI.J, U.CI.J, Length.J)

c(L.CI.HPD.J,  U.CI.HPD.J, Length.HPD.J)

c(L.CI.U,  U.CI.U, Length.U)

c(L.CI.HPD.U, U.CI.HPD.U, Length.HPD.U)

c(L.CI.F, U.CI.F, Length.F)

##############################################################################