# 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[0]
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, 2021 at 16:12
• Either you are unlucky or you're doing something wrong. A guess: samples.std is incorrect.
– whuber
Jul 28, 2021 at 16:16
• @whuber ahhhh yeah.... there's an $n-1$ supposed to be in there somewhere!
– Him
Jul 28, 2021 at 16:18
• @whuber (0.9485, 0.9508) NICE!
– Him
Jul 28, 2021 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, 2021 at 13:51
• @MattKrause the robot tells me: "You can accept your own answer tomorrow"
– Him
Jul 29, 2021 at 15:08
• Opps, forgot about the time-out! Just didn't want you to be bashful about accepting your own answer. Jul 29, 2021 at 23:33

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)
mean(CI[1,] <= 0 & CI[2,] >= 0)
[1]  0.94907
sd(CI[1,]<=0 & CI[2,]>=0)/sqrt(10^5)
[1] 0.0006952454  # aprx 95% margin of simulation error

• Wow, this is so..... concise. :D I see you did this as 100,000 bernoulli trials. I did mine as 1000 n=1000 binomial trials. This way makes a lot more sense, and probably results in the tightest CI on the coverage prob....
– Him
Jul 28, 2021 at 20:52
• 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, 2021 at 21:36