13

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))
milo
  • 103
  • 3
Him
  • 2,027
  • 10
  • 25

2 Answers2

19

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)

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

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
BruceET
  • 47,896
  • 2
  • 28
  • 76
  • 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 '21 at 20:52
  • 4
    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. – BruceET Jul 28 '21 at 21:36