# 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...?

Thank you :)

-
A quick note related to your question about coverage probability: We have some unknown quantity that we cannot observe, in this case, the mean. We can construct a confidence interval, which gives us a rough idea of where the true value of the parameter is. Depending on what points we sample from in our population, we will get different confidence intervals. Coverage probability is the proportion of those intervals that contain the true value of the parameter. See Wikipedia: en.wikipedia.org/wiki/Coverage_probability – Christopher Aden Mar 22 '11 at 2:32

## 2 Answers

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.

-
Thank you so much for this - it makes sense :). I get ~0.92 as my answer. – user3832 Mar 22 '11 at 11:00

You have several issues with your code:

1. 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.
2. (since solved) You have UCL less than LCL, which is a slightly odd use of upper and lower
3. You do not need semicolons in R unless you want more than one instruction on the same line
4. (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.
5. 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)

-
Ok - I goofed on the UCL and LCL values. I'm used to java code with semi-colons. I fixed the post but still get 0... :( Could you please elaborate on what you mean by the mean of the empty set? Sorry for my bad english. – user3832 Mar 22 '11 at 0:50
Instead of telling me everything that wrong with current code, could you give hint at how to approach from beginning? – user3832 Mar 22 '11 at 1:31
The zeros stem from the fact that mean(LCL > 0 & UCL < 0); should be mean(LCL < 2 & UCL > 2) (the number of times the Ci includes 2, the true mean, divided by m). – JMS Mar 22 '11 at 2:43