##############################################################################
# R code to apply with rainfall data 
# For the Difference Mean of Delta-gamma Distributions (Equal 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 2

D2 =	c(0,56.8,66,0,0,0.8,4,1.6,0,0
		,0,1,0,2.2,2.9,0,0,0,0,0
		,0,0,0,13.6,0,1.9,0,47.3,16.5,47.3
		,0,0,24.4,0,6.2,9.2,33.9,21.5,0,60.1
		,32.3,47.3,1,71.5,0,0)

  x 	<- D2^(1/3)
  n_x 	<- length(D2)
  
  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 3

D3 =	c(3.8,0,38.3,0,0,0,0,0,4.7,0
		,0,19,12.2,0,21.8,0,23.5,0,0,0
		,31.3,0,0,38.8,32.4,0,4.1,0,0.7,0
		,0,0,8.3,0,0,0,0.5,15.3,0,0
		,0.6,0,0,0,0,50.2)

  y 	<- D3^(1/3)
  n_y 	<- length(D3)
  
  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)

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