根据样本数据计算置信区间

Posted 2023-03-12

【中文标题】根据样本数据计算置信区间

【英文标题】：Compute a confidence interval from sample data

【发布时间】：2013-02-08 14:55:14

【问题描述】：

我有一些样本数据，我想为其计算置信区间，假设为正态分布。

我已经找到并安装了 numpy 和 scipy 包，并且已经让 numpy 返回平均值和标准差（numpy.mean(data)，数据为列表）。任何有关获取样本置信区间的建议将不胜感激。

【问题讨论】：

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

对于非分类数据需要做什么，例如回归任意实数值？

回答我自己的问题，是的，请参阅：stats.stackexchange.com/questions/554332/…

【参考方案1】：

关于 Ulrich 的回答 - 即使用 t 值。当真正的方差未知时，我们使用它。这是您唯一拥有的数据是示例数据的情况。

对于 bogatron 的回答，这涉及到 z 表。在已知并提供方差时使用 z 表。然后你也有样本数据。 Sigma 不是样本均值的估计标准差。已经知道了。

假设方差是已知的，我们想要 95% 的置信度：

from scipy.stats import norm
alpha = 0.95
# Define our z
ci = alpha + (1-alpha)/2
#Lower Interval, where n is sample siz
c_lb = sample_mean - norm.ppf(ci)*((sigma/(n**0.5)))
c_ub = sample_mean + norm.ppf(ci)*((sigma/(n**0.5)))

只有样本数据和未知方差（这意味着必须仅根据样本数据计算方差），Ulrich 的答案非常有效。但是，您可能希望指定置信区间。如果您的数据是 a 并且您希望置信区间为 0.95：

import statsmodels.stats.api as sms
conf = sms.DescrStatsW(a).tconfint_mean(alpha=0.05)
conf

【参考方案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 或Wikipedia）。

对于大样本量 n，样本均值呈正态分布，可以使用 st.norm.interval() 计算其置信区间（如 Jaime 的评论中所建议的那样）。但上述解决方案对于小 n 也是正确的，其中st.norm.interval() 给出的置信区间太窄（即“假置信”）。有关更多详细信息，请参阅我的answer 到类似问题（以及此处的 Russ 的 cmets 之一）。

这是一个示例，其中正确的选项给出（基本上）相同的置信区间：

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)

【讨论】：

我相信您应该致电 st.t.interval(0.05) 以获得 95% 的置信区间。

不，st.t.interval(0.95) 对于 95% 置信区间是正确的，请参阅 docs 以获得 scipy.stats.t。不过，SciPy 将参数命名为 alpha 似乎不太理想。

如果我有两个数据数组，然后计算它们的平均值之差。有什么方法可以得到这个平均差的 95% CI？您能想出任何简单的方法来像您在此处提供的那样使用 StatsModelsl 吗？

当样本量较小（小于30）时，应使用Student-t分布，在本例中为（[10,11,12,13）。结果，正态分布给出了不同的结果。例如，如果您将样本量增加到 1000，则 t- 和 norm 给出几乎相同的结果。

这对分类和回归有用吗？ （例如，如果有负值，任意大小）

【参考方案3】：

从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.from_samples(data)) 创建一个 NormalDist 对象，这使我们可以通过 NormalDist.mean 和 NormalDist.stdev 访问样本的均值和标准差。

使用累积分布函数的倒数 (inv_cdf)，根据给定置信度的标准正态分布（由 NormalDist() 表示）计算 Z-score。

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

---

这假设样本量足够大（假设超过 ~100 个点），以便使用标准正态分布而不是学生的 t 分布来计算 z 值。

【讨论】：

这对分类和回归有用吗？ （例如，如果有负值，任意大小）

是否有任何理由使用错误但近似正确的正态分布而不是完全正确的 t 分布？我看不出使用正确的 t 分布（见***.com/a/15034143/7735095 或***.com/a/34474255/7735095）与观察次数无关。

【参考方案4】：

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

你可以这样计算。

【讨论】：

小心“私人”使用sp.stats.t._ppf。如果没有进一步的解释，我对此并不满意。最好直接使用sp.stats.t.ppf，除非您确定自己知道自己在做什么。在快速检查the source 时，_ppf 跳过了相当多的代码。可能是良性的，但也可能是不安全的优化尝试？

这对分类和回归有用吗？ （例如，如果有负值，任意大小）

回答自己：是的，因为它正在计算 CI mu_n +- t.val(0.95) * std_n/ sqrt(n) 详情请参阅：stats.stackexchange.com/questions/554332/…

【参考方案5】：

首先从look-up table 中查找z-value 以获得所需的置信区间。置信区间为mean +/- z*sigma，其中sigma 是样本均值的估计标准偏差，由sigma = s / sqrt(n) 给出，其中s 是根据样本数据计算的标准偏差，n 是样本大小.

【讨论】：

scipy.stats.norm.interval(confidence, loc=mean, scale=sigma)

@bogatron，关于置信区间的建议微积分，不会是 mean +/- z * sigma/sqrt(n)，其中 n 是样本大小？

@David，你是对的。我弄错了sigma 的含义。我的回答中的sigma 应该是样本均值的估计标准差，而不是分布的估计标准差。我已经更新了答案以澄清这一点。感谢您指出这一点。

@Jaime 评论中有误导性。如果您正在计算 t 学生置信区间，则不使用 sigma，而是使用标准误差 sigma/np.sqrt（观察总数），否则会得到错误的结果。你也可以说：scipy.stats.norm.interval(confidence, loc=mean, scale=standard error)

“查看查找表”是此堆栈交换的不适当答案。它应该是库调用的一部分，以便代码可以在运行时自行获取 z 分数，并且可以将置信区间作为变量公开给用户。