从样本数据计算置信区间

Question

我有一些样本数据，想要计算一个置信区间，假设这些数据符合正态分布。

我已经找到了并安装了numpy和scipy这两个包，并且用numpy计算出了平均值和标准差（用numpy.mean(data)，其中data是一个列表）。如果有人能给我一些关于如何获取样本置信区间的建议，我将非常感激。

Answer 1

从 Python 3.8 开始，标准库里提供了一个叫 NormalDist 的对象，它是 statistics 模块的一部分：

from statistics import NormalDist

def confidence_interval(data, confidence=0.95):
  dist = NormalDist.from_samples(data)
  z = NormalDist().inv_cdf((1 + confidence) / 2.)
  h = dist.stdev * z / ((len(data) - 1) ** .5)
  return dist.mean - h, dist.mean + h

这个对象：

• 可以从数据样本中创建一个 NormalDist 对象（使用 NormalDist.from_samples(data)），这样我们就能通过 NormalDist.mean 和 NormalDist.stdev 来获取样本的平均值和标准差。

• 可以根据标准正态分布（用 NormalDist() 表示）计算给定置信度的 Z-score，这通过累积分布函数的反函数来实现（使用 inv_cdf）。

• 根据样本的标准差和平均值，生成置信区间。

---

这里假设样本量足够大（比如超过100个数据点），这样才能使用标准正态分布，而不是学生的t分布来计算 z 值。

Answer 2

这里是shasan代码的简化版本，用来计算数组a的平均值的95%置信区间：

import numpy as np, scipy.stats as st

st.t.interval(0.95, len(a)-1, loc=np.mean(a), scale=st.sem(a))

不过，使用StatsModels的tconfint_mean方法可能更好：

import statsmodels.stats.api as sms

sms.DescrStatsW(a).tconfint_mean()

这两种方法的基本假设是，样本（数组a）是从一个标准差未知的正态分布中独立抽取的（可以参考MathWorld或维基百科）。

对于样本量较大的情况，样本的平均值是正态分布的，可以使用st.norm.interval()来计算它的置信区间（正如Jaime的评论中提到的）。但是，上述方法在样本量小的情况下也是正确的，因为st.norm.interval()会给出过于狭窄的置信区间（也就是“假信心”）。有关更多细节，可以查看我对类似问题的回答（以及Russ在这里的评论）。

这里有一个例子，显示正确的方法给出的置信区间（基本上）是相同的：

In [9]: a = range(10,14)

In [10]: mean_confidence_interval(a)
Out[10]: (11.5, 9.4457397432391215, 13.554260256760879)

In [11]: st.t.interval(0.95, len(a)-1, loc=np.mean(a), scale=st.sem(a))
Out[11]: (9.4457397432391215, 13.554260256760879)

In [12]: sms.DescrStatsW(a).tconfint_mean()
Out[12]: (9.4457397432391197, 13.55426025676088)

最后，使用st.norm.interval()得到的错误结果：

In [13]: st.norm.interval(0.95, loc=np.mean(a), scale=st.sem(a))
Out[13]: (10.23484868811834, 12.76515131188166)

Answer 3

import numpy as np
import scipy.stats

def mean_confidence_interval(data, confidence=0.95):
    a = 1.0 * np.array(data)
    n = len(a)
    m, se = np.mean(a), scipy.stats.sem(a)
    h = se * scipy.stats.t.ppf((1 + confidence) / 2., n-1)
    return m, m-h, m+h

你可以这样计算。