python - 用 scipy 获取置信区间的正确方法

Question

我有一个一维数据数组:

a = np.array([1,2,3,4,4,4,5,5,5,5,4,4,4,6,7,8])

我想获得 68% 的置信区间(即:1 sigma)。

this answer 中的第一条评论声明这可以使用 scipy.stats.norm 中的 scipy.stats.norm.interval 来实现函数，通过:

from scipy import stats
import numpy as np
mean, sigma = np.mean(a), np.std(a)

conf_int = stats.norm.interval(0.68, loc=mean, 
    scale=sigma)

但是 this post 中的评论指出获得置信区间的实际正确方法是:

conf_int = stats.norm.interval(0.68, loc=mean, 
    scale=sigma / np.sqrt(len(a)))

即sigma除以样本量的平方根:np.sqrt(len(a))。

问题是:哪个版本是正确的？

Answer 1

单次抽签的 68% 置信区间来自正态分布 平均 mu 和标准偏差 sigma 是

stats.norm.interval(0.68, loc=mu, scale=sigma)

N 的平均值的 68% 置信区间来自正态分布 平均 mu 和标准偏差 sigma 是

stats.norm.interval(0.68, loc=mu, scale=sigma/sqrt(N))

直观地说，这些公式是有道理的，因为如果你拿着一 jar 软糖 bean 让很多人猜测软糖 bean 的数量，每个人可能会相差很多——相同的标准偏差 sigma -- 但是猜测的平均值在估计实际数字方面会做得非常好，这反射(reflect)在均值的标准偏差缩小了 1/sqrt(N ).

---

如果单次平局有方差sigma**2，则由Bienaymé formula ，N uncorrelated 绘制的总和具有方差 N*sigma**2。

均值等于总和除以 N。当您将随机变量(如总和)乘以常数时，方差乘以常数的平方。那是

Var(cX) = c**2 * Var(X)

所以均值的方差等于

(variance of the sum)/N**2 = N * sigma**2 / N**2 = sigma**2 / N

因此均值的标准差(即方差的平方根)等于

sigma/sqrt(N).

这是分母中sqrt(N)的原点。

---

以下是一些示例代码，基于 Tom 的代码，用于演示上述声明:

import numpy as np
from scipy import stats

N = 10000
a = np.random.normal(0, 1, N)
mean, sigma = a.mean(), a.std(ddof=1)
conf_int_a = stats.norm.interval(0.68, loc=mean, scale=sigma)

print('{:0.2%} of the single draws are in conf_int_a'
      .format(((a >= conf_int_a[0]) & (a < conf_int_a[1])).sum() / float(N)))

M = 1000
b = np.random.normal(0, 1, (N, M)).mean(axis=1)
conf_int_b = stats.norm.interval(0.68, loc=0, scale=1 / np.sqrt(M))
print('{:0.2%} of the means are in conf_int_b'
      .format(((b >= conf_int_b[0]) & (b < conf_int_b[1])).sum() / float(N)))

打印

68.03% of the single draws are in conf_int_a
67.78% of the means are in conf_int_b

请注意，如果您使用 mean 和 sigma 的估计值定义 conf_int_b根据示例 a，平均值可能不会落在 conf_int_b 中频率。

---

如果您从分布中获取一个样本并计算样本均值和标准差，

mean, sigma = a.mean(), a.std()

请注意，不能保证这些将等于 人口 均值和标准差，并且我们假设人口是正态分布的——这些不是自动给定的!

如果您抽样并想要估计总体均值和标准偏差，你应该使用

mean, sigma = a.mean(), a.std(ddof=1)

因为 sigma 的这个值是 unbiased estimator为总体标准差。