如何用 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 设计用于拟合给定样本的正态分布参数(mu 和 sigma) : NormalFactory()
创建此 class.
的一个实例
方法 build
进行实际拟合,returns class Normal
的对象表示参数为 mu[ 的正态分布=44=] 和 sigma 从样本估计。
sample
整形是为了确保 OpenTURNS 理解输入 sample
是 one-dimension 点的集合,而不是单个 multi-dimensional 点。
然后 class Normal
提供方法 computeQuantile
来计算分布的任何分位数(本例中为第 95 个百分位数)。
此解决方案不计算精确的公差界限,因为它使用正态分布的分位数而不是学生 t-distribution。实际上,这意味着它忽略了 mu 和 sigma 上的估计误差。实际上,这只是样本量非常小的问题。
为了说明这一点,这里比较了标准正态 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
我想根据已知 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
.
import openturns as ot
ot.NormalFactory().build(sample.reshape(-1, 1)).computeQuantile(0.95)
让我们打开包装。 NormalFactory
是一个 class 设计用于拟合给定样本的正态分布参数(mu 和 sigma) : NormalFactory()
创建此 class.
方法 build
进行实际拟合,returns class Normal
的对象表示参数为 mu[ 的正态分布=44=] 和 sigma 从样本估计。
sample
整形是为了确保 OpenTURNS 理解输入 sample
是 one-dimension 点的集合,而不是单个 multi-dimensional 点。
然后 class Normal
提供方法 computeQuantile
来计算分布的任何分位数(本例中为第 95 个百分位数)。
此解决方案不计算精确的公差界限,因为它使用正态分布的分位数而不是学生 t-distribution。实际上,这意味着它忽略了 mu 和 sigma 上的估计误差。实际上,这只是样本量非常小的问题。
为了说明这一点,这里比较了标准正态 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