############################################################################## # R code to apply with rainfall data # For the Single Mean of Delta-gamma Distributions ############################################################################## library(stats) library(invgamma) library(EnvStats) library(HDInterval) library(extraDistr) library(MASS) library(fitdistrplus) library(LaplacesDemon) # Single # x1 = data # M = number of simulation runs # m = number of generalized computation # n0 = size of zero in population # n1 = size of non-zero in population # shape = shape parameter # rate = rate parameter # tau = mean # 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 <- rep(0,m) ex.jef <- rep(0,m) var.jef <- rep(0,m) delta.jef <- rep(0,m) tau.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 <- rep(0,m) ex.uni <- rep(0,m) var.uni <- rep(0,m) delta.uni <- rep(0,m) tau.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 <- rep(0,m) beta1 <- rep(0,m) R.delta <- rep(0,m) R.mu <- rep(0,m) R.sigsq <- rep(0,m) R.mean <- rep(0,m) R.tau <- rep(0,m) L.CI.F <- numeric() U.CI.F <- numeric() CP.F <- numeric() Length.F <- numeric() y_MLE <- matrix(rep(0,2)) ############################################################################## # Start Iteration ############################################################################## #Dataset 1 D1 = c(0,0,0,74.4,0,63.3,36.7,0,0,0 ,0,0,0,10.9,0,0,0,0,4.8,0 ,0,0.3,0,0,0,52.9,7.2,30.5,0,0 ,3.6,20.9,1.9,0,0,0,16.4,0,12,2.8 ,10.4,15.1,0,82,31,29.6,0.9,41.8,0,21.7) x <- D1^(1/3) n <- length(D1) nonzero_value <- x[which(x!=0)] n_nonzero <- length(nonzero_value) n_zero <- n-n_nonzero y <- nonzero_value y_MLE <- fitdistr(y, "gamma", start=list(shape=1, rate=1))$estimate shape <- y_MLE[1] rate <- y_MLE[2] scale <- 1/rate delta <- n_zero/(n_zero+n_nonzero) tau <- (1-delta)*shape/rate alpha <- 0.05 mu_x <- (shape/scale)^1/3*(1-(1/(9*shape))) sigsq_x <- 1/(9*(shape)^(1/3)*(scale^(2/3))) xbar <- mean(nonzero_value) s <- sd(nonzero_value) sum_y <- sum((y-mean(y))^2) ex_hat <- xbar var_hat <- s^2 u <- rchisq(m , n_nonzero-1)/(n_nonzero-1) z <- rnorm(m) mu <- xbar + z*s/sqrt(n_nonzero)/sqrt(u) sigsq <- s^2/u ############################################################################## # compute CI based on credible confidence interval Jeffreys prior ############################################################################## beta.jef <- rbeta(m,(n-n_nonzero)+(1/2),(n_nonzero)+(3/2)) ex.jef <- rnorm(m,ex_hat,sigsq_x/n_nonzero) var.jef <- rinvgamma(m, n_nonzero/2,(sum_y)/2) delta.jef <- beta.jef mean.jef <- (0.5*ex.jef + sqrt(0.25*(ex.jef)^2 + var.jef))^3 tau.jef <- delta.jef*mean.jef L.CI.J <- quantile(tau.jef,alpha/2) U.CI.J <- quantile(tau.jef,(1-alpha/2)) Length.J <- U.CI.J- L.CI.J ############################################################################## # compute CI based on HPD Jeffreys prior ############################################################################## HPD.J <- hdi(tau.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 <- rbeta(m,(n-n_nonzero)+1,(n_nonzero)+1) ex.uni <- rnorm(m,ex_hat,sigsq_x/n_nonzero) var.uni <- rinvgamma(m, (n_nonzero-3)/2,(sum_y)/2) delta.uni <- beta.uni mean.uni <- (0.5*ex.uni + sqrt(0.25*(ex.uni)^2 + var.uni))^3 tau.uni <- delta.uni*mean.uni L.CI.U <- quantile(tau.uni,alpha/2) U.CI.U <- quantile(tau.uni,(1-alpha/2)) Length.U <- U.CI.U - L.CI.U ############################################################################## # compute CI based on HPD uniform prior ############################################################################## HPD.U <- hdi(tau.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 <- rbeta(m,n_nonzero,n_zero+1) beta1 <- rbeta(m,n_nonzero+1,n_zero) R.delta <- (1/2)*beta0 + (1/2)*beta1 R.mu <- xbar + z*s/sqrt(n_nonzero)/sqrt(u) R.sigsq <- s^2/u R.mean <- (0.5*R.mu + sqrt(0.25*(R.mu)^2 + R.sigsq))^3 R.tau <- R.delta*R.mean L.CI.F <- quantile(R.tau,alpha/2,type=8,na.rm = TRUE) U.CI.F <- quantile(R.tau,(1-alpha/2),type=8,na.rm = TRUE) Length.F <- U.CI.F - L.CI.F ############################################################################## c(n, delta, shape, rate, tau) 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) ##############################################################################