我正在模拟一组从正态分布中提取的样本的 95% 置信区间。既然数据是正常的，那么我认为，我 95% 的置信度应该转化为 95% 的覆盖概率。但是，我得到了大约 94% 的结果。具体来说，我采用 1000 个大小为 n=10 的样本来制作一堆 CI 并估计覆盖概率，然后重复1000次以获得覆盖概率的 CI。我的覆盖概率的 5 sigma CI 原来是 ~(0.9384, 0.9408)。是否有一些统计原因，或者我做错了什么？

这是我的模拟代码：

    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))