# Not getting 95% coverage for 95% t-distribution CI

I'm simulating a bunch of 95% confidence intervals on samples taken from a normal distribution. Since the data is normal, then, I think, my 95% confidence should translate into a 95% coverage probability. However, I'm getting something like 94%. Specifically, I'm taking 1000 samples of size n=10 to make a bunch of CIs and estimate a coverage probability, then doing that 1000 times to get a CI for the coverage probability. My five sigma CI for the coverage probability turns out to be ~(0.9384, 0.9408). Is there some statistical reason for this, or am I doing something wrong?

Here's my simulation code:

    import numpy as np
import scipy.stats as stats

def CI_coverage(alpha, dist, n, n_samples):
''' creates n_samples samples of size n
creates an 1-alpha confidence interval for each
computes the fraction of those that contain mu '''
# get samples
samples = np.stack([dist.rvs(size=n) for i in range(n_samples)])

# summary stats
mu = dist.mean()
xbar = samples.mean(axis=1)
s = samples.std(axis=1)

# compute CIs... note that xbar, s, CI_low, CI_high are arrays size n_samples
t = stats.t.ppf(1 - alpha/2, n-1)
interval_width = t * s / np.sqrt(n)
CI_low = xbar - interval_width
CI_high = xbar + interval_width

coverage_p = np.sum(np.logical_and(CI_low < mu, mu < CI_high)) / samples.shape
return coverage_p

mu = 1
sigma = 0.5
norm_dist = stats.norm(loc=mu, scale=sigma)

n = 10
n_samples = 1000
n_CI_samples = 1000
# compute the empirical coverage probability many times
CI_coverages = [CI_coverage(0.05, norm_dist, n, n_samples) for i in range(n_CI_samples)]

# use this to get a CI for the coverage probabilities
CI_c_mean = np.mean(CI_coverages)
CI_c_std = np.std(CI_coverages)

print(CI_c_mean - 5*CI_c_std / np.sqrt(n_CI_samples), CI_c_mean + 5*CI_c_std / np.sqrt(n_CI_samples))

• I didn't check your Python code for glitches. However, 1000 iterations is not enough to tell the difference between .94 and .95; try 100,000. Jul 28 at 16:12
• Either you are unlucky or you're doing something wrong. A guess: samples.std is incorrect.
– whuber
Jul 28 at 16:16
• @whuber ahhhh yeah.... there's an $n-1$ supposed to be in there somewhere!
– Him
Jul 28 at 16:18
• @whuber (0.9485, 0.9508) NICE!
– Him
Jul 28 at 16:20

Per @whuber's comment, np.std() provides a biased estimate of the sample standard deviation. Fortunately, the function allows you to correct for that by specifying a number of degrees of freedom with the ddof parameter:

s = samples.std(axis=1, ddof=1)


Fixing this gives coverage probabilities that are consistent with the expected 95% CI: (0.9485, 0.9508)

• Glad you figured it out! You can, if you want, "accept" your own answer. You won't get any reputation from it, but it'll indicate to future readers that the question is "resolved". Jul 29 at 13:51
• @MattKrause the robot tells me: "You can accept your own answer tomorrow"
– Him
Jul 29 at 15:08
In R, using $-notation to pick the 95% CI out of t.test output, I get $$0.949 \pm 0.001,$$ from 100,000 iterations. set.seed(2021) CI = replicate( 10^5, t.test(rnorm(10))$conf.int)

• Using a built-in procedure in R such as t.test makes the code shorter and ensures you are using a well-vetted method. However, running time can be somwhat longer than for more direct code because R does all the work of formatting complete output, of which only a relatively small part is used. Jul 28 at 21:36