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

直觉上，这些公式是有道理的，因为如果你举起一罐果冻豆，让很多人猜猜果冻豆的数量，每一个个体都可能有很大的偏离——同样的标准差sigma——但是猜测的平均值在估计实际数量方面会做得非常好，这反映在平均收缩的标准差乘以系数1/sqrt(N)上。

---

如果一次抽签有方差sigma**2，那么通过Bienaymé formula，不相关的N抽签的和有方差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。

Answer 2

我用一个已知置信区间的数组测试了你的方法。normal（mu，std，size）返回一个以mu为中心的数组，标准偏差为std（在the docs中，定义为Standard deviation (spread or “width”) of the distribution.）。

from scipy import stats
import numpy as np
from numpy import random
a = random.normal(0,1,10000)
mean, sigma = np.mean(a), np.std(a)
conf_int_a = stats.norm.interval(0.68, loc=mean, scale=sigma)
conf_int_b = stats.norm.interval(0.68, loc=mean, scale=sigma / np.sqrt(len(a)))


conf_int_a
(-1.0011149125527312, 1.0059797764202412)
conf_int_b
(-0.0076030415111100983, 0.012467905378619625)

由于sigma值应为-1到1，因此/ np.sqrt(len(a))方法似乎不正确。

编辑

由于我没有资格发表上述评论，我将澄清这一回答如何与联合国大学的彻底答复联系起来。如果用正态分布填充一个随机数组，总数的68%将落在平均值的1-σ内。在上面的例子中，如果你检查一下

b = a[np.where((a>-1)&(a <1))]
len(a)
> 6781

或者68%的人口在1σ之内。嗯，大约68%。当你使用一个越来越大的数组时，你将接近68%（在10次试验中，9次在-1和1之间）。这是因为1-σ是数据的固有分布，数据越多，就越能解决它。

基本上，我对你问题的解释是如果我有一个数据样本，我想用它来描述它们的分布，那么找到这些数据的标准差的方法是什么？而联合国大学的解释似乎更为我可以用68%的置信度将平均值置于什么区间？。这意味着，对于果冻豆，我回答了他们是如何猜测的，而unutbu回答了他们的猜测告诉我们关于果冻豆的什么。

Answer 3

我刚刚检查了R和GraphPad是如何计算置信区间的，如果样本量（n）很小，它们会增加置信区间。E、 g，对于n=2，比一个大n要大6倍以上。这个代码（基于shasan的answer）与它们的置信区间相匹配：

import numpy as np, scipy.stats as st

# returns confidence interval of mean
def confIntMean(a, conf=0.95):
  mean, sem, m = np.mean(a), st.sem(a), st.t.ppf((1+conf)/2., len(a)-1)
  return mean - m*sem, mean + m*sem

对于R，我检查了t检验（a）。GraphPad的confidence interval of a mean页面有关于样本大小依赖关系的“用户级”信息。

这里是Gabriel示例的输出：

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

In [3]: confIntMean(a, 0.68)
Out[3]: (3.9974214366806184, 4.877578563319382)

In [4]: st.norm.interval(0.68, loc=np.mean(a), scale=st.sem(a))
Out[4]: (4.0120010966037407, 4.8629989033962593)

注意，confIntMean()和st.norm.interval()间隔之间的差异在这里相对较小；len（a）==16不是太小。