################################################################################# # R code to construct confidence intervals # ################################################################################# library(HDInterval) library(invgamma) CI.CVLN <- function(M,m,k,n1,n2,n3,mean,sqrt.var1,sqrt.var2,sqrt.var3) { # begin process # x = data # ni = number of sample size # M = number of simulation runs # m = number of generalized computation # s = sample standard deviation # sigma = population standard deviation # theta = coefficient of variation # alpha = level of significance # L.CI1 = lower CI based on FGCI # U.CI1 = upper CI based on FGCI # L.CI2 = lower CI based on MOVER # U.CI2 = upper CI based on MOVER # L.CI3 = lower CI based on CA # U.CI3 = upper CI based on CA # L.CI4 = lower CI based on Bayesian # U.CI4 = upper CI based on Bayesian # Ex.length1 = average length of CI based on FGCI # Ex.length2 = average length of CI based on MOVER # Ex.length3 = average length of CI based on CA # Ex.length4 = average length of CI based on Bayesian # CP.FGCI = coverage probability of CI based on GCI # CP.MOVER = coverage probability of CI based on MOVER # CP.CA = coverage probability of CI based on CA # CP.BS = coverage probability of CI based on Bayesian ni <- rep(0,k) s <- rep(0,k) xbar <- rep(0,k) Rtheta.FGCI <- rep(0,m) mu.hat.RML <- rep(0,k) sigma.hat.RML <- rep(0,k) theta.CA <- rep(0,m) theta.BS <- rep(0,m) CP.FGCI <- rep(0,M) CP.MOVER <- rep(0,M) CP.CA <- rep(0,M) CP.BS <- rep(0,M) Length.FGCI <- rep(0,M) Length.MOVER <- rep(0,M) Length.CA <- rep(0,M) Length.BS <- rep(0,M) alpha <- 0.05 # fix level of significance z.alpha <- qnorm(1-(alpha/2)) # compute z at 1-alpha/2 sqrt.var <- c(sqrt.var1,sqrt.var2,sqrt.var3) theta1 <- sqrt(exp(sqrt.var1^2)-1) # compute population CV theta2 <- sqrt(exp(sqrt.var2^2)-1) # compute population CV theta3 <- sqrt(exp(sqrt.var3^2)-1) # compute population CV var.theta1 <- ((sqrt.var1^4)*(exp(2*(sqrt.var1^2))))/(2*(n1-1)*((exp(sqrt.var1^2))-1)) var.theta2 <- ((sqrt.var2^4)*(exp(2*(sqrt.var2^2))))/(2*(n2-1)*((exp(sqrt.var2^2))-1)) var.theta3 <- ((sqrt.var3^4)*(exp(2*(sqrt.var3^2))))/(2*(n3-1)*((exp(sqrt.var3^2))-1)) theta.p <- c(theta1,theta2,theta3) var.theta.p <- c(var.theta1,var.theta2,var.theta3) theta <- (sum(theta.p/var.theta.p))/(sum(1/var.theta.p)) #------------------------------------------------------------------------------------------------------------------- # Start iteration #------------------------------------------------------------------------------------------------------------------- for (i in 1:M) { # begin loop i #------------------------------------------------------------------------------------------------------------------- # Generate x #------------------------------------------------------------------------------------------------------------------- x1 <- rnorm(n1,mean,sqrt.var1) # generate x from normal xbar.1 <- mean(x1) # compute sample mean of xi s.1 <- sd(x1) # compute standard deviation of xi x2 <- rnorm(n2,mean,sqrt.var2) # generate x from normal xbar.2 <- mean(x2) # compute sample mean of xi s.2 <- sd(x2) # compute standard deviation of xi x3 <- rnorm(n3,mean,sqrt.var3) # generate x from normal xbar.3 <- mean(x3) # compute sample mean of xi s.3 <- sd(x3) # compute standard deviation of xi #---------------------------------------------------------------------------------------------------------------- # Call data #----------------------------------------------------------------------------------------------------------------- ni <- c(n1,n2,n3) # call sample size xbar <- c(xbar.1,xbar.2,xbar.3) # call mean s <- c(s.1,s.2,s.3) # call standard deviation (sample) thetahat <- sqrt((exp(s^2))-1) # compute CV using sample variance var.thetahat <- ((s^4)*(exp(2*(s^2))))/(2*(ni-1)*((exp(s^2))-1)) # compute variance of CV frac1 <- sum(thetahat/var.thetahat) frac2 <- sum(1/var.thetahat) thetahat.large <- frac1/frac2 # compute CV of LN #------------------------------------------------------------------------------------------------------------------- # Compute CI based on FGCI #------------------------------------------------------------------------------------------------------------------- for (j in 1:m) { # begin loop j x1.star <- rnorm(n1,mean,sqrt.var1) # generate x1.star from normal distribution xbar1.star <- mean(x1.star) # compute sample mean of x1.star s1.star <- sd(x1.star) # compute standard deviation of x1.star x2.star <- rnorm(n2,mean,sqrt.var2) # generate x2.star from normal distribution xbar2.star <- mean(x2.star) # compute sample mean of x2.star s2.star <- sd(x2.star) # compute standard deviation of x2.star x3.star <- rnorm(n3,mean,sqrt.var3) # generate x3.star from normal distribution xbar3.star <- mean(x3.star) # compute sample mean of x3.star s3.star <- sd(x3.star) # compute standard deviation of x3.star xbar.star <- c(xbar1.star,xbar2.star,xbar3.star) # call mean s.star <- c(s1.star,s2.star,s3.star) # call standard deviation (sample) Rsig.sqrt.FGCI <- (s^2)*(sqrt.var^2)/(s.star^2) R.FGCI <- sqrt(exp(Rsig.sqrt.FGCI)-1) # compute FGPQ for CV of LN Rvar.FGCI <- ((Rsig.sqrt.FGCI^2)*(exp(2*Rsig.sqrt.FGCI)))/(2*(ni-1)*((exp(Rsig.sqrt.FGCI))-1)) Rtheta.FGCI[j] <- sum(R.FGCI/Rvar.FGCI)/sum(1/Rvar.FGCI) } # end loop j L.CI1 <- quantile(Rtheta.FGCI,0.025,type=8) # compute lower bound U.CI1 <- quantile(Rtheta.FGCI,0.975,type=8) # compute upper bound CP.FGCI[i] <- ifelse(L.CI1 0) stop("Root is not bracketed in the specified interval \n") chg <- upper-lower while (abs(chg) > tol) { theta.bisect.new <- (lower+upper)/2 f.new <- fn(theta.bisect.new,ni,s) if (abs(f.new) <= tol) break if (f.lo*f.new < 0) upper <- theta.bisect.new if (f.hi*f.new < 0) lower <- theta.bisect.new chg <- upper-lower feval <- feval+1 } list(theta.bisect.new,chg,feval) } data1 <-bisect(fn1, 0, 2, theta.bisect =theta.bisect, ni=ni, s=s) theta.hat.RML <- data1[[1]] mu.hat.RML <- xbar mu.hat.RML1 <- mu.hat.RML[[1]] mu.hat.RML2 <- mu.hat.RML[[2]] mu.hat.RML3 <- mu.hat.RML[[3]] sigma.hat.RML <- s sigma.hat.RML1 <- sigma.hat.RML[[1]] sigma.hat.RML2 <- sigma.hat.RML[[2]] sigma.hat.RML3 <- sigma.hat.RML[[3]] #----------------------------------------------------------------------------------------------------------------- # Generate artificial sample #----------------------------------------------------------------------------------------------------------------- for (j in 1:m) { # begin loop j x.RML1 <- rnorm(n1,mu.hat.RML1,sigma.hat.RML1) # generate x from normal xbar.RML.1 <- mean(x.RML1) # compute sample mean of xi s.RML.1 <- sd(x.RML1) # compute standard deviation of xi x.RML2 <- rnorm(n2,mu.hat.RML2,sigma.hat.RML2) # generate x from normal xbar.RML.2 <- mean(x.RML2) # compute sample mean of xi s.RML.2 <- sd(x.RML2) # compute standard deviation of xi x.RML3 <- rnorm(n3,mu.hat.RML3,sigma.hat.RML3) # generate x from normal xbar.RML.3 <- mean(x.RML3) # compute sample mean of xi s.RML.3 <- sd(x.RML3) # compute standard deviation of xi #----------------------------------------------------------------------------------------------------------------- # Call data #----------------------------------------------------------------------------------------------------------------- xbar.RML <- c(xbar.RML.1,xbar.RML.2,xbar.RML.3) # call mean s.RML1 <- c(s.RML.1,s.RML.2,s.RML.3) # call standard deviation thetahat.RML <- sqrt((exp(s.RML1^2))-1) var.RML <- ((s.RML1^4)*(exp(2*(s.RML1^2))))/(2*(ni-1)*((exp(s.RML1^2))-1)) frac1.RML <- sum(thetahat.RML/var.RML) frac2.RML <- sum(1/var.RML) theta.CA[j] <- frac1.RML/ frac2.RML # compute CV of LN } # end loop j L.CI3 <- quantile(theta.CA,0.025,type=8) # compute lower bound U.CI3 <- quantile(theta.CA,0.975,type=8) # compute upper bound CP.CA[i] <- ifelse(L.CI3