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))
  • 2
    $\begingroup$ 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. $\endgroup$
    – BruceET
    Jul 28 at 16:12
  • 6
    $\begingroup$ Either you are unlucky or you're doing something wrong. A guess: samples.std is incorrect. $\endgroup$
    – whuber
    Jul 28 at 16:16
  • 2
    $\begingroup$ @whuber ahhhh yeah.... there's an $n-1$ supposed to be in there somewhere! $\endgroup$
    – Him
    Jul 28 at 16:18
  • $\begingroup$ @whuber (0.9485, 0.9508) NICE! $\endgroup$
    – 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)

  • 3
    $\begingroup$ 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". $\endgroup$ Jul 29 at 13:51
  • $\begingroup$ @MattKrause the robot tells me: "You can accept your own answer tomorrow" $\endgroup$
    – Him
    Jul 29 at 15:08
  • 1
    $\begingroup$ Opps, forgot about the time-out! Just didn't want you to be bashful about accepting your own answer. $\endgroup$ Jul 29 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.

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
  • $\begingroup$ 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.... $\endgroup$
    – Him
    Jul 28 at 20:52
  • 4
    $\begingroup$ 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. $\endgroup$
    – BruceET
    Jul 28 at 21:36

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.