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?
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
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 <= mean <= interval: true += 1 print true/nsim