# Monte Carlo comparison of performances of different confidence intervals # for a distribution mean --- exponential distribution setting. set.seed(321) # I'll initialize the results matrix. results <- matrix( rep( 0, 9 ), nrow=3 ) rownames(results) <- c("Student's t","percentile","bootstrap t") colnames(results) <- c("miss left","miss right","miss total") results # I'll go through a loop 2500 times. On each pass through the loop I'll # generate a new data vector, compute 3 c.i's, and for each one determine # it it misses the estimand by falling to the left of it or falling to the # right of it. for (i in 1:2500) { x <- rexp(25, 1) t.result <- t.test( x, conf.level=0.90 ) if ( t.result$conf.int[1] > 1 ) { results[1,2] <- results[1,2] + 1 } else { if ( t.result$conf.int[2] < 1 ) { results[1,1] <- results[1,1] + 1 } } bss <- matrix( sample( x, 25*999, replace=T ), nrow=999 ) mean.rep <- apply( bss, 1, mean ) sorted.mean.rep <- sort( mean.rep ) if ( sorted.mean.rep[50] > 1 ) { results[2,2] <- results[2,2] + 1 } else { if ( sorted.mean.rep[950] < 1 ) { results[2,1] <- results[2,1] + 1 } } sd.rep <- apply( bss, 1, sd ) mean.x <- mean( x ) z.rep <- ( mean.rep - mean.x )/( sd.rep/5 ) sorted.z.rep <- sort( z.rep ) se.hat <- sd( x )/5 if ( mean.x - sorted.z.rep[950]*se.hat > 1 ) { results[3,2] <- results[3,2] + 1 } else { if ( mean.x - sorted.z.rep[50]*se.hat < 1 ) { results[3,1] <- results[3,1] + 1 } } } results[,3] <- results[,1] + results[,2] results/2500