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

Posted 2023-03-12

【中文标题】用 scipy 获取置信区间的正确方法

【英文标题】：Correct way to obtain confidence interval with scipy

【发布时间】：2015-03-30 07:47:12

【问题描述】：

我有一个一维数据数组：

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

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

【问题讨论】：

您可以使用引导程序：***.com/a/66008548/10375049

【参考方案1】：

我使用具有已知置信区间的数组测试了您的方法。 numpy.random.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)) 方法似乎不正确。

编辑

由于我没有在上面发表评论的声誉，我将澄清这个答案如何与 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 回答 他们的猜测告诉了我们关于果冻豆的什么信息。

【参考方案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))

直观地说，这些公式是有道理的，因为如果你举起一罐果冻豆并让很多人猜测果冻豆的数量，每个人可能会相差很多——相同的标准偏差@ 987654325@ -- 但猜测的平均值在估计实际数字方面会做得非常好，这反映在平均值的标准差缩小了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。

【讨论】：

很好的答案@unutbu，非常彻底。只是为了确保我已经涵盖了我的基础，你如何解释汤姆在另一个答案中给出的例子？

我认为汤姆和我的回答是一致的。我根据他的显示添加了一些代码，表明单次抽签和均值确实以预期的频率落在置信区间内。

scipy.stats 有 sem 函数 stats.sem(a) == a.std(ddof=1) / np.sqrt(len(a)) 除了浮点错误

R 和 Graphpad 在计算置信区间时使用与样本量相关的调整，这对于小样本量有显着差异。详情见我的回答。

如果单个样本已被引导 N 次，是否还应该使用 scale=sigma/sqrt(N)？

【参考方案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 页面有关于样本大小依赖性的“用户级别”信息。

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