如何用 scipy 计算单侧公差区间

How to calculate one-sided tolerance interval with scipy

我想根据已知 N(样本大小)、标准差和平均值的数据集的正态分布计算单侧公差界限。

如果间隔是双向的,我会执行以下操作:

conf_int = stats.norm.interval(alpha, loc=mean, scale=sigma)

在我的情况下,我正在引导示例,但如果不是,我会在 Whosebug 上参考此 post: 并使用以下内容:conf_int = stats.norm.interval(0.68, loc=mean, scale=sigma / np.sqrt(len(a)))

你会如何做同样的事情,但要将其计算为单边边界(95% 的值高于或低于 x<--bound)?

我假设您有兴趣使用正态分布计算 one-side 公差界限(基于您提到 scipy.stats.norm.interval 函数作为您需要的 two-sided 等价物这一事实) .

那么好消息是,基于 tolerance interval Wikipedia page:

One-sided normal tolerance intervals have an exact solution in terms of the sample mean and sample variance based on the noncentral t-distribution.

(仅供参考:不幸的是,two-sided 设置并非如此)

此断言基于this paper。此外第 4.8 段(第 23 页)提供了公式。

坏消息是,我认为没有 ready-to-use scipy 函数可以安全地调整并用于您的目的。

但是你可以很容易地自己计算出来。您可以在 Github 存储库中找到包含此类计算器的存储库,您可以从中找到灵感,例如 that one 我从中构建了以下说明性示例:

import numpy as np
from scipy.stats import norm, nct

# sample size
n=1000

# Percentile for the TI to estimate
p=0.9
# confidence level
g = 0.95

# a demo sample
x = np.array([np.random.normal(100) for k in range(n)])

# mean estimate based on the sample
mu_est = x.mean()

# standard deviation estimated based on the sample
sigma_est = x.std(ddof=1)

# (100*p)th percentile of the standard normal distribution
zp = norm.ppf(p)

# gth quantile of a non-central t distribution
# with n-1 degrees of freedom and non-centrality parameter np.sqrt(n)*zp
t = nct.ppf(g, df=n-1., nc=np.sqrt(n)*zp)

# k factor from Young et al paper
k = t / np.sqrt(n)

# One-sided tolerance upper bound
conf_upper_bound = mu_est + (k*sigma_est)

这是一个使用 openturns 库的 one-line 解决方案,假设您的数据是一个名为 sample.

的 numpy 数组
import openturns as ot
ot.NormalFactory().build(sample.reshape(-1, 1)).computeQuantile(0.95)

让我们打开包装。 NormalFactory 是一个 class 设计用于拟合给定样本的正态分布参数(musigma) : NormalFactory() 创建此 class.

的一个实例

方法 build 进行实际拟合,returns class Normal 的对象表示参数为 mu[ 的正态分布=44=] 和 sigma 从样本估计。

sample 整形是为了确保 OpenTURNS 理解输入 sample 是 one-dimension 点的集合,而不是单个 multi-dimensional 点。

然后 class Normal 提供方法 computeQuantile 来计算分布的任何分位数(本例中为第 95 个百分位数)。

此解决方案不计算精确的公差界限,因为它使用正态分布的分位数而不是学生 t-distribution。实际上,这意味着它忽略了 musigma 上的估计误差。实际上,这只是样本量非常小的问题。

为了说明这一点,这里比较了标准正态 N(0,1) 分布的 PDF 和具有 19 个自由度的学生 t-distribution 的 PDF(这意味着样本量为20).他们几乎无法区分。

deg_freedom = 19
graph = ot.Normal().drawPDF()
student = ot.Student(deg_freedom).drawPDF().getDrawable(0)
student.setColor('blue')
graph.add(student)
graph.setLegends(['Normal(0,1)', 't-dist k={}'.format(deg_freedom)])
graph