基于最小样本空间的Johnson分布拟合方法

张 琪 吴义忠 刘 鑫 乔 平

1.华中科技大学国家数控系统工程技术研究中心，武汉，4300742.华中科技大学国家企业信息化应用支撑软件工程技术研究中心，武汉，430074

0 引言

现代复杂机电产品设计常受到知识的缺乏和产品运行环境的变化等方面不确定性因素影响，致使设计的机电产品在运行时部分可靠度指标可能会发生变化或偏移导致引发产品故障[1]，因此，在产品设计优化阶段需基于大量的样本数据来充分考虑各方面的不确定性因素。然而，实际优化过程常受到计算资源的制约,在满足特定条件下应尽可能地减少估值采样的次数。根据少量样本点拟合出总体分布，基于该分布计算可靠度指标，可有效避免多次昂贵估值带来的巨大计算损耗，提高计算效率。

在对总体分布未知的样本数据进行分布函数拟合时，可选择自适应的分布形式对总体进行拟合，如标准的二参数和三参数Weibull分布，调整其参数即可实现分布模型接近于指数分布、正态分布等分布模型。在样本数据较多的情况下，图解法、极大似然估计法[2]、最小二乘法关于水平残差和垂直残差平方和最小[3]等方法对所需参数的估值较为精确。近年来，Weibull分布在机电产品的可靠性分析中有着非常重要的应用[2-5]，然而Weibull分布中只有两个或三个参数，不能全面表达总体的分布特征，同时，根据有限的样本数据获得所需参数的精确估值还没有得到有效的解决[6]。

四参数的Johnson分布族函数[7](下文简称Johnson分布)于1949年被提出，后来发展成包含四种不同类型的分布族函数。由于具有多参数、多类型的特征，Johnson分布通过对样本拟合可表达出更多的总体的分布特征。Johnson分布具有较强的自适应性，选择合适的类型和调整参数，可接近于任意一个标准连续型分布模型，其中包括Weibull分布模型[8]。DEBROTA等[9]提出配矩法(moment matching)、百分位数配比法(percentile matching)、最小二乘法(least squares)和最小范数估值法(minimum Lpnorm estimation)四种方法，这些方法根据有限的样本数据即可实现对Johnson分布中四个参数的精确估计。

本文引用Johnson分布作为未知总体的分布形式对样本进行拟合，并基于假设检验[10]，控制第Ⅰ类错误和第Ⅱ类错误的发生概率，根据Z检验法中双边检验问题的施行特征函数和曲线，推导出样本空间的最小值。最后设计实验进行一致性检验，采用K-S(Kolmogorov-Smirnov)拟合检验法检验总体的拟合分布与真实分布的一致性。实验结果表明，当样本空间达到最小样本空间时，Johnson分布根据有限样本拟合出的总体分布与样本的真实分布具有一致性，满足拟合精度要求。

1 Johnson分布族函数及其拟合

1.1 Johnson分布族函数

为便于统一表达各种不同类型的连续型随机变量的累积分布函数FX(x)=Pr(X≤x)和概率密度函数f=F′(x)，Johnson提出了四参数的Johnson分布族函数。Johnson分布的概率密度函数为

(1)

其中，γ和δ为形状参数；ε为位置参数；λ为尺度参数；f(·)为简单的函数表达式，根据f(·)的不同，Johnson分布可分为对数正态(lognormal, SL)、无界 (unbounded, SU)、有界(bounded, SB)和正态(normal, SN)四种不同类型。

考虑到客观条件的制约，生产实践中获得的数据大多数情况下都服从SB类型的Johnson分布。以SB类型为例，三类参数对Johnson分布模型的影响如图1～图3所示。

图1 形状参数γ、δ对模型的影响Fig.1 The influence of shape parameters γ、δon the model

图2 位置参数ε对模型的影响Fig.2 The influence of location parameters εon the model

图3 尺度参数λ对模型的影响Fig.3 The influence of dimension parameters λon the model

f(·)和f′(·)的具体表达式如下：

x的取值范围H为

对式(1)进行积分，即可得到如下各个类型相应的累积分布函数：

无界型累积分布函数

-∞

正态型累积分布函数

-∞

对数正态型累积分布函数

FSL(x)=

有界型累积分布函数

FSB(x)=

式中，Erf(·)为高斯误差函数；Erfc(·)为误差互补函数；ArcSinh(·)为反双曲正弦函数。

Johnson分布通过变形可得

(2)

该式为Johnson转换式，可将连续型随机变量X映射到一个服从标准正态分布的随机变量Z上[7]。通过Johnson转换式可将样本数据{x1,x2,…,xn}转换成一组服从标准正态分布的新样本{z1,z2,…,zn}，根据样本对总体X进行Johnson分布拟合时，拟合精度越高, 则新样本{z1,z2,…,zn}对标准正态分布的服从度就越高。

1.2 Johnson分布拟合

对一组样本的总体分布进行拟合时，首先根据样本所反映的总体分布特征选择合适的拟合分布形式，再估算出分布形式中各个参数的值，从而得到总体的拟合分布。而Johnson分布具有较强的自适应性，通过选择类型和调整参数，Johnson分布模型可近似于其他不同形式的连续型分布模型，即Johnson分布可统一表达不同形式的连续型分布。因此，当对一组真实分布形式完全未知的样本进行拟合时，不需要分析样本所反映的总体分布特征来选择合适的分布形式进行总体拟合，直接选用Johnson分布作为总体的分布形式加以拟合即可。

本文通过MATLAB软件编程实现了Johnson分布对总体拟合的过程。根据样本的分布特征，基于Johnson分布类型的选择准则正确选择出合适的分布类型，并得到其相应的参数函数表达式，即概率密度函数表达式和累积分布函数表达式。考虑到百分位数配比法比其他方法对所需参数进行估值时更加简便且能保证估值精度[11]，本文采用百分位数配比法对所需参数进行精确估计，将各个未知参数的估值代入参数函数表达式中即可得到总体的拟合分布。其拟合流程如图4所示。

图4 Johnson分布对总体拟合流程图Fig.4 The flow chart of Johnson distribution fitting

2 最小样本空间分析

基于假设检验，通过控制第Ⅰ类错误和第Ⅱ类错误发生的概率，根据Z检验法中双边检验问题的施行特征函数和曲线，在保证拟合精度的条件下对拟合所需的最小样本空间进行推导。

基于Johnson分布对随机变量X的总体分布进行拟合后得到的是复杂的分布函数表达式，若直接根据该分布进行假设检验推导出拟合所需最小的样本空间则更为复杂。文献[12]提出一种简便的解决方法：通过Fisher转换将随机变量映射到一个服从正态分布的随机变量上，在正态分布上推导出拟合所需的最小样本空间。由于对总体进行Johnson分布拟合的精度与通过转换式(2)得到的新样本对标准正态分布的服从度相关，所以本文采用Johnson转换将随机变量X映射到一个服从标准正态分布的随机变量Z上，在标准正态分布上通过控制第Ⅰ类错误和第Ⅱ类错误发生的概率得到样本空间的范围。为此，引入施行特征函数[10]。

定义若C是参数θ的某检验问题的一个检验法,则称β(θ)=Pθ(接受H0)为检验法C的施行特征函数或OC函数，其图形称为OC曲线。

根据上述定义，考虑到双边检验问题H0：μ=μ0,H1：μ≠μ0, 正态总体均值的Z检验法的OC函数为

(3)

式中，Φ(·)为求标准正态分布的累计概率密度函数。

其OC曲线如图5所示，β(μ)是|λ|的严格单调下降函数。

图5 OC曲线图Fig.5 OC curve graph

在双边检验问题中，若要求对H1中满足|μ-μ0|≥δ>0的μ处的函数值β(μ)≤β，则需要解超越方程：

由此知只要样本空间n满足

即只要n满足

(4)

就能使当μ∈H1且|μ-μ0|≥δ(δ>0，为取定的值)时，第Ⅰ类错误发生的概率不超过给定的值α，第Ⅱ类错误发生的概率不超过给定的值β。

当通过式(4)确定样本空间范围时，由文献[13]可知，容许误差δ是假设检验所试图揭示样本与整体之间的差异大小，总体标准差σ体现个体变异度，对于没有给定专业意义上的容许误差水平的情况，用0.25倍或0.50倍的σ来设定δ。本文取α=0.05，β=0.05，δ=0.5σ。查表得zα/2=z0.025=1.96，zβ=z0.05=1.645，将上述数据代入式(4)得

即通过控制第Ⅰ类错误和第Ⅱ类错误发生的概率可得到最小的样本空间为52。

3 实验验证与结果分析

基于假设检验推导出当样本空间不小于52时，采用Johnson分布对总体的分布进行拟合可达到理论精度。考虑到样本空间较小且变量为连续型时，K-S拟合检验法比χ2拟合检验法具有更强的检出力，因此本文采用K-S检验法设计实验来验证当样本空间为52时，Johnson分布对总体的拟合是否能够达到要求的精度。实验步骤如下：

(1)采集样本点。从真实分布已知的总体X中随机采集4组样本, 各组样本空间分别为42，47，52，57。

(2)确定分布的表达式。将步骤(1)作为Johnson拟合程序的输入，运行程序, 输出即为步骤(1)中各组样本所对应的Johnson分布,记拟合分布的总体为X。当X～N(0,1)、X～Γ(1,1)、X～W(1,1)且样本空间为52时，Johnson分布拟合效果如图6～图8所示。

图6 真实总体分布为X～N(0,1)Fig.6 The real population distribution subjectin g to X～N(0,1)

图7 真实总体分布为X～Γ(2,2)Fig.7 The real population distribution subjectin g to X～Γ(2,2)

图8 真实总体分布为X～W(3,5)Fig.8 The real population distribution subjectin g to X～W(3,5)

图9 X～N时，K-S检验结果Fig.9 The result of K-S test, if X～N

图10 X～Γ时，K-S检验结果Fig.10 The result of K-S test, if X～Γ

图11 X～W时，K-S检验结果Fig.11 The result of K-S test, if X～W

本实验中，步骤(1)中总体的真实分布形式选取为几种常见的连续型分布，包括正态分布、Gamma分布、Weibull分布等, 实验结果如图9～图 11所示。根据图9～图11所示的实验检验结果，通过分析可知:

(1)不失一般性，随着样本空间增大，拟合精度提高，K-S检验法的检验结果q值呈增大趋势。当样本空间为52时，检验结果q值均大于0.05，说明假设成立，从X和X中分别随机抽取的两组样本来自于同一分布。这表明，样本空间为52时，总体的拟合分布与真实分布具有一致性，即在一定精度条件下，拟合分布等价于真实分布。

(2)当样本空间为52，样本数据取自真实分布X～N(2,2)和X～W(3,5)时，K-S检验法得到的q值远大于0.05，说明实际拟合精度远高于要求精度。但一般情况下，q值均以0.05为下界，在一定范围内波动。这表明，采用Johnson分布对一组总体分布形式完全未知的样本进行拟合时，52作为最小样本空间具有普适性。

4 工程算例

曲轴是发动机中最重要的部件，它的可靠性对发动机的使用寿命有着重要的影响。当基于应力-强度干涉模型计算曲柄颈的可靠度指标时，需在工作过程中统计曲柄颈危险截面处应力值的变化从而获得其概率密度。将曲柄臂抽象为矩形，如图12所示，当曲轴各参数取表1中的值时，已知曲柄颈处所受的力F服从正态分布N(15,1)，即可根据上述数据构建曲轴的力学模型，进而求得所需应力值。

图12 曲轴力学模型Fig.12 The mechanical model of crankshaft

F(kN)N(15,1)W(N·m)6.6l1(mm)400l2(mm)220l3(mm)100e(mm)120d(mm)60φ(°)13

利用测力传感器，从服从正态分布N(15,1)的F中随机抽取52个值作为样本数据点，通过力学分析和公式计算[14]，曲柄颈危险截面处对应的52个应力值如表2所示。

采用Johnson分布对所求的52个曲柄颈处所受力F进行拟合，得到应力F的概率密度函数

采用Johnson分布对所求的52个应力值进行拟合，得到应力σ0的概率密度函数为

表2 52个F样本点及其对应的危险截面处应力σ0

力F和应力σ0的概率分布拟合如图13和图14所示。

图13 力F的概率分布拟合Fig.13 The fitting pdf of F

图14 应力σ0的概率分布拟合Fig.14 The fitting pdf of σ0

5 结论

通过引入Johnson分布函数，采用百分位数配比法求解Johnson分布函数的参数，并基于假设检验原理确定满足精度要求的最小样本空间，同时通过K-S拟合检验法对总体真实分布与拟合分布进行一致性验证，保证拟合精度满足要求，实现了Johnson分布对总体分布拟合过程。本文的研究工作可快速确定对总体分布拟合所需的最小样本空间，提高基于少量昂贵估值对总体分布拟合的精度和效率，为可靠性设计优化过程中的可靠度计算提供理论依据，对提高可靠性设计优化的计算效率具有重要意义。