使用 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 states that this can be achieved using scipy.stats.norm.interval from the scipy.stats.norm函数中的第一条评论，来自：

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

我使用已知置信区间的数组测试了您的方法。 numpy.random.normal(mu,std,size) returns 一个以 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)) 方法似乎不正确。

编辑

由于我没有资格在上面发表评论，所以我将阐明这个答案如何与 unutbu 的详尽答案联系起来。如果您使用正态分布填充随机数组，则总数的 68% 将落在均值的 1-σ 范围内。在上面的例子中，如果你检查你看到

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

或 68% 的人口落在 1σ 范围内。嗯，大约 68%。随着您使用越来越大的阵列，您将接近 68%（在 10 个试验中，9 个在 -1 和 1 之间）。那是因为 1-σ 是数据的固有分布，你拥有的数据越多，你就越能解决它。

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

Answer 2

一次抽取的正态分布的 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))

直觉上，这些公式是有道理的，因为如果你举起一罐糖豆并让很多人猜测糖豆的数量，每个人可能会相差很多——同样的标准deviation 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 3

我刚刚检查了 R 和 GraphPad 如何计算置信区间，它们在样本量 (n) 较小的情况下增加了区间。例如，与大 n 相比，n=2 的结果超过 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.test(a)。 GraphPad 的 confidence interval of a mean 页面有 "user level" 关于样本量依赖性的信息。

这里是 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 不算小了。

使用 scipy 获得置信区间的正确方法

Correct way to obtain confidence interval with scipy

python

numpy

scipy

confidence-interval

编辑