n <- 25         # set n to desired sample size 
M <- 1000000    # set M to desired number of Monte Carlo trials
set.seed(13)    # setting seed for random number generator, so can reproduce later (using the same seed)

df <- 2*n 
ql <- qchisq(0.05, df)
qu <- qchisq(0.95, df) 
z <- qnorm(0.95) 
zn <- z/sqrt(n) 

c1 <- 0 
c2 <- 0  
c3 <- 0 
c4 <- 0 
c5 <- 0 
w1 <- numeric(M) 
w2 <- numeric(M) 
w3 <- numeric(M) 
w4 <- numeric(M) 

for (m in 1:M) { 
  x <- rexp(n, 1)    # generate sample of n observations from exponential dist'n having a mean of 1 
  
  xbar <- mean(x) 
  exact <- 2*n*xbar 
  l1 <- exact/qu 
  u1 <- exact/ql 
  if ( (l1 < 1) && (u1 > 1) ) c1 <- c1 + 1 
  w1[m] <- u1 - l1
  
  l2 <- xbar*(1 - zn) 
  u2 <- xbar*(1 + zn) 
  if ( (l2 < 1) && (u2 > 1) ) c2 <- c2 + 1
  w2[m] <- u2 - l2 

  l3 <- xbar/(1 + zn) 
  u3 <- xbar/(1 - zn) 
  if ( (l3 < 1) && (u3 > 1) ) c3 <- c3 + 1
  w3[m] <- u3 - l3 

  exppart <- exp(zn) 
  l4 <- xbar/exppart 
  u4 <- xbar*exppart 
  if ( (l4 < 1) && (u4 > 1) ) c4 <- c4 + 1 
  w4[m] <- u4 - l4 

  if ( ((l3 + l4)/2 < 1) && ((u3 +u4)/2 > 1) ) c5 <- c5 + 1 
}
 
# exact interval, est. cov. prob. & average width 
c1/M 
mean(w1) 

# 1st approx. interval, est. cov. prob. & average width
c2/M 
mean(w2) 

# 2nd approx. interval, est. cov. prob. & average width
c3/M 
mean(w3) 

# 3rd approx. interval, est. cov. prob. & average width 
c4/M 
mean(w4) 

# averaging conf. bounds for 2nd and 3rd approx. estimators, cov. prob. 
c5/M