# Can we be really confident about confidence intervals?

I created a population of a normally distributed value. Then I pulled 100 000 samples out of it. For each of them I calculated 95% confidence z-interval. I expected to see that at least 95% of intervals would contain a population mean, but in fact only about 94% did. I tried this many times, and there was 94% again and again. I also tried t-intervals with the same result. How can it be explained?

• The method you are using doesn't actually achieve 95% confidence. Use the correct method, which employs the Student t distribution. You will definitely not get the same result: if you do, there must be an error in the calculation. – whuber Feb 24 '17 at 21:20
• @whuber You were right, I've tried t-intervals from Scipy and they are correct. But why that first method doesn't achieve 95% confidence? – Michael Boojum Feb 24 '17 at 21:46
• I am unaware of any authoritative reference that claims it does give a 95% confidence interval. It's only an approximation that was superseded when William Gosset proposed using the Student t distribution in a 1908 paper. – whuber Feb 24 '17 at 22:09
• @whuber Thank you! Your answers are great as always. – Michael Boojum Feb 24 '17 at 22:11
• Just a note: in practice, a method that claims 95% confidence but achieves 94% would be, in my opinion, great. This is because all models have a assumptions required to achieve this 95% level. The question is not whether your data meets these assumptions, but how badly it fails these assumptions. If these assumptions are only mildly inaccurate (or, moreover, breaking these assumptions only mildly affects your inference), maybe your method results in 93% true confidence level. If these assumptions are terrible, maybe you have 63% true confidence level. – Cliff AB Feb 24 '17 at 23:04

The zconfint() function takes its critical values from the standard normal distribution, instead of from the $t$ distribution (which is appropriate for small samples such as yours). The normal critical values tend to be smaller, which yields narrower confidence intervals, and explains why your intervals seem to succeed in covering less frequently. If you increase the size argument from 30 to something like 300 then you should see coverage closer to (but still slightly less than) 0.95.

If you really want to construct a $t$ confidence interval, you should create a DescrStatsW object and call the tconfint_mean() function:

import numpy as np
from statsmodels.stats.weightstats import DescrStatsW

pop = np.random.randn(1000000)

true = 0.
mean = pop.mean()
nsim = 10000
for i in range(nsim):
sample = np.random.choice(pop, size=30)
d = DescrStatsW(sample)
interval = d.tconfint_mean()
if interval[0] <= mean <= interval[1]:
true += 1

print true/nsim