Take the 2-minute tour ×
Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. It's 100% free, no registration required.

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 :)

share|improve this question
    
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 2

up vote 4 down vote accepted

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.

share|improve this answer
    
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)
share|improve this answer
    
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

Your Answer

 
discard

By posting your answer, you agree to the privacy policy and terms of service.