Monte Carlo experiment to estimate coverage probability I'm working on a problem as follows for a course that I'm auditing:

Suppose a 95% symmetric t-interval is applied to estimate a mean, but the
  sample data are non-normal. Then the probability that the confidence interval covers
  the mean is not necessarily equal to 0.95.
  Use a Monte Carlo experiment to
  estimate the coverage probability of
  the t-interval for random samples of
  $\chi^2(2)$ data with sample size $n = 20$.

Here is the current state of my R code:
alpha = 0.05;
n = 20;
m = 1000;

UCL = numeric(m);
LCL = numeric(m);

for(i in 1:m)
{
    x = rchisq(n, 2);
    LCL[i] = mean(x) - qt(alpha / 2, lower.tail = FALSE) * sd(x);
    UCL[i] = mean(x) + qt(alpha / 2, lower.tail = FALSE) * sd(x);
}

# This line below is wrong...
mean(LCL > 0 & UCL < 0);

The problem is that the result is $0$. Am I approaching this question incorrectly? What exactly does coverage probability mean...?
 A: I disagree with Henry - I think you should be dividing by sqrt(n), because it's a confidence interval for the mean. You also have to add a df = n-1 argument to your qt calls.
And the last line should be mean(LCL < 2 & UCL > 2). This is because 2 is the true mean, and you're interested in the condition that 2 is in the confidence interval.
A: You have several issues with your code:


*

*Your mean(UCL < 0 & LCL > 0) is decidedly strange, and in particular is failing because UCL is coming out positive so you are taking the mean of an empty set.  A $\chi^2$ distribution takes only positive values.

*(since solved) You have UCL less than LCL, which is a slightly odd use of upper and lower

*You do not need semicolons in R unless you want more than one instruction on the same line 

*(Wrong - as pointed out by mark999) You are dividing by sqrt(n).  This wrongly narrows your confidence intervals: it is for finding the standard error of the mean, but you care about the original distribution.

*The question tells you to use "t-interval" but you are using a normal distribution. You might try typing ?qt into R


Try this
alpha = 0.05
n = 20
m = 1000

UCL = numeric(m)
LCL = numeric(m)

for(i in 1:m)
{
    x = rchisq(n, 2) # compare with x = rnorm(n) + 2
    LCL[i] = mean(x) - qt(alpha / 2, df=n-1, lower.tail = FALSE)*sd(x)/sqrt(n)
    UCL[i] = mean(x) + qt(alpha / 2, df=n-1, lower.tail = FALSE)*sd(x)/sqrt(n)
}

mean(LCL < 2 & UCL > 2)

