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

Question

我有一个1维数据数组:

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

我希望获得68%置信区间(即:1西格玛).

在第一个评论这个回答指出,这可以实现使用scipy.stats.norm.interval从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)

但是这篇文章中的评论指出,获得置信区间的实际正确方法是:

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

也就是说,sigma除以样本大小的平方根:np.sqrt(len(a)).

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

Answer 1

具有平均μ和std偏差sigma的正态分布的单次抽取的68%置信区间是

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

N的平均值的68%置信区间来自正态分布,平均μ和std偏差sigma为

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

直觉上,这些公式是有道理的,因为如果你拿起一罐果冻豆并要求大量的人猜测果冻豆的数量,每个人可能会大量减少 - 相同的标准偏差sigma- 但是估计的平均值将非常精确地估计实际数量,这可以通过平均收缩的标准差来反映1/sqrt(N).

如果单个绘制具有方差sigma**2,则通过Bienaymé公式,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

请注意,如果您conf_int_b使用估计值mean并sigma 基于样本进行定义a,则均值可能不会conf_int_b与所需频率一致.

如果你把样品从分布和计算样本均值和标准方差,

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

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

如果您采样并想要估算总体均值和标准差,则应使用

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

因为西格玛的这个值是人口标准差的无偏估计.

Answer 2

我使用具有已知置信区间的数组测试了您的方法.numpy.random.normal(mu,std,size)返回一个以mu为中心的数组,标准偏差为std(在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 3

我刚刚检查了R和GraphPad如何计算置信区间,并且在样本量较小的情况下它们会增加间隔(n).例如,与大n相比,n = 2超过6倍.此代码(基于shasan的答案)匹配其置信区间:

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的平均页面的置信区间具有关于样本大小依赖性的"用户级别"信息.

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