如何使用scipy计算单边公差区间

问题描述

在给定具有已知N（样本大小），标准差和均值的数据集的情况下，我想根据正态分布计算一个单边公差范围。

如果间隔是双面的，我将执行以下操作：

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

在我的情况下，我正在引导示例，但如果不是，我将在stackoverflow上参考这篇帖子：Correct way to obtain confidence interval with scipy并使用以下内容：conf_int = stats.norm.interval(0.68,scale=sigma / np.sqrt(len(a)))

您将如何做同样的事情，但是将其计算为单边绑定（95％的值高于或低于x

解决方法

我假设您对使用正态分布计算单边公差范围感兴趣（基于您提到scipy.stats.norm.interval函数是您需要的两面等效的事实）。

那么，好消息是，基于tolerance interval Wikipedia page：

单侧法线公差区间在样本均值和基于非中心t分布的样本方差方面都有精确的解决方案。

（仅供参考：不幸的是，双面设置不是这种情况）

此断言基于this paper。除了第4.8节（第23页）之外，还提供了公式。

坏消息是，我认为没有可以安全调整和使用的现成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库的单行解决方案，假设您的数据是一个名为sample的numpy数组。

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

让我们打开包装。 NormalFactory是一个类，旨在适合给定样本的正态分布（ mu 和 sigma ）的参数：NormalFactory()创建一个这个班。

方法build进行实际拟合，并返回类Normal的对象，该对象表示参数 mu 和 sigma 的正态分布从样本中估算出来。

在那里进行sample重塑，以确保OpenTURNS理解输入sample是一个一维点的集合，而不是一个多维点。

类Normal然后提供方法computeQuantile来计算分布的任何分位数（在此示例中为第95个百分位数）。

此解决方案未计算确切的公差范围，因为它使用了来自正态分布的分位数而不是Student t分布。实际上，这意味着它忽略了 mu 和 sigma 的估计误差。实际上，这仅是很小样本量的问题。

为说明这一点，这里是标准正态N（0,1）分布的PDF与自由度为19的Student t分布的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