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

张 琪1 吴义忠2 刘 鑫2 乔 平2

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

摘要基于应力-强度干涉模型计算可靠度指标时,常使用统计方法获得应力和强度的概率分布,然而,在复杂机电产品的设计过程中,每一个样本值的获取都是一次昂贵估值,大量样本值的获取面临耗时严重等问题。针对该问题,采用参数化的Johnson分布拟合概率分布,研究基于最小样本空间的Johnson分布拟合方法。采用百分位数配比法求解Johnson分布族函数的参数,并基于假设检验原理确定满足拟合精度要求的最小样本空间,同时通过K-S拟合检验法对总体真实分布与拟合分布进行一致性分析,保证拟合精度符合要求。工程应用结果表明,该方法基于少量的昂贵估值样本可快速得到样本的总体分布,为可靠度指标的快速计算奠定了基础。

关键词Johnson分布族函数;假设检验;最小样本空间;一致性分析

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 Lp norm estimation)四种方法,这些方法根据有限的样本数据即可实现对Johnson分布中四个参数的精确估计。

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

1 Johnson分布族函数及其拟合

1.1 Johnson分布族函数

为便于统一表达各种不同类型的连续型随机变量的累积分布函数FX(x)=Pr(Xx)和概率密度函数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)进行积分,即可得到如下各个类型相应的累积分布函数:

无界型累积分布函数

-∞<x<+∞

正态型累积分布函数

-∞<x<+∞

对数正态型累积分布函数

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函数为

β(μ)=Pμ(接受H0)=


Φ(zα/2-λ)-Φ(-zα/2-λ)=
Φ(zα/2-λ)+Φ(zα/2+λ)-1

(3)

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

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

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

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

解出上述方程可确定样本容量n。通常因n值较大,故总可以认为于是故有

由此知只要样本空间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。当XN(0,1)、XΓ(1,1)、XW(1,1)且样本空间为52时,Johnson分布拟合效果如图6~图8所示。

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

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

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

(3)一致性检验。从真实分布X和拟合分布X中分别随机抽取一组样本空间为1 000的样本,假设这两组样本来自于同一分布,并采用K-S双样本检验法对该假设加以检验。实现检验过程中选择临界α=0.05,计算样本累计频率与理论分布累计概率的绝对差,令最大的绝对差为与临界值作比较并返回到q值。当q>α时,假设成立,这两组样本来自于同一分布,即样本的真实分布X和拟合分布X具有一致性。反之假设不成立。

图9 XN时,K-S检验结果
Fig.9 The result of K-S test, if XN

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

图11 XW时,K-S检验结果
Fig.11 The result of K-S test, if XW

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

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

(2)当样本空间为52,样本数据取自真实分布XN(2,2)和XW(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

表1 曲轴参数

Tab.1 The parameters 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

Tab.2 52 samples of F and the stress σ0 at th e corresponding critical section

序号F(kN)截面应力σ0(MPa)114.956 18110.385 3214.714 35108.726 2315.825 22116.354 1413.666 32101.546 3514.168 63104.985 3614.225 53105.375 2714.727 53108.816 6813.422 9499.881 6914.455 47106.951 01014.111 97104.597 21115.303 52112.769 81214.979 97110.548 61314.160 41104.929 01414.251 12105.550 51513.843 60102.759 51614.646 15108.258 41714.997 15110.666 51813.792 15102.407 41915.350 18113.090 22016.354 59119.994 42115.353 02113.109 72214.020 79103.972 72315.701 54115.504 02414.466 44107.026 22515.964 23117.309 72615.337 56113.003 52715.489 97114.050 32815.520 06114.257 12915.840 38116.458 23015.269 54112.536 43116.378 97120.162 13214.519 35107.388 93313.039 1097.258 43414.399 67106.568 53515.739 36115.764 03613.941 82103.431 93714.406 36106.614 43814.722 13108.779 53914.531 38107.471 54014.802 30109.329 54113.927 84103.336 24213.249 7998.697 94314.700 93108.634 24415.507 97114.174 04514.721 94108.778 24615.022 89110.843 14714.738 00108.888 44815.025 55110.861 44915.281 98112.621 95015.033 48110.915 85114.422 91106.727 85216.127 49118.432 3

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分布对总体分布拟合过程。本文的研究工作可快速确定对总体分布拟合所需的最小样本空间,提高基于少量昂贵估值对总体分布拟合的精度和效率,为可靠性设计优化过程中的可靠度计算提供理论依据,对提高可靠性设计优化的计算效率具有重要意义。

参考文献

[1] 陈小前,姚雯,欧阳琦. 飞行器不确定性多学科设计优化理论与应用[M].北京: 科学出版社,2013:1-10.

CHEN Xiaoqian, YAO Wen, OUYANG Qi. Theory and Application of Uncertainty-based Multidisciplinary Design Optimization for Flight Vehicles[M]. Beijing: Science Press, 2013:1-10.

[2] XIE M, TANG Y, GOH T N. A Modified Weibull Extension with Bathtub-shaped Failure Rate Function[J]. Reliability Engineering & System Safety, 2002, 76(3): 279-285.

[3] 丛伟,陈晓阳,王志坚,等. Weibull分布产品小样本定时截尾试验方案下的可靠性评估[J].中国机械工程,2013,24(14):1891-1896

CONG Wei, CHEN Xiaoyang, WANG Zhijian, et al. Reliability Evaluation of Products with Life Following Weibull Distribution under Type I Censoring Test Plan with Small Sample Size[J]. China Mechanical Engineering, 2013, 24(14): 1891-1896.

[4] 高顺川,冯静,孙权,等.基于威布尔分布的动态故障树定量分析方法[J]. 质量与可靠性,2005(5):32-35.

GAO Shunchuan, FENG Jing, SUN Quan, et al. Dynamic Fault Tree Quantitative Analysis Method Based on Weibull Distribution[J]. Quality and Reliability, 2005(5):32-35.

[5] 郑锐. 三参数威布尔分布参数估计及在可靠性分析中的应用[J]. 振动与冲击, 2015, 34(5):78-81.

ZHENG Rui. Parameter Estimation of Three Parameter Weibull Distribution and Its Application in Reliability Analysis[J]. Journal of Vibration and Shock, 2015, 34(5):78-81.

[6] 凌丹. 威布尔分布模型及其在机械可靠性中的应用研究[D]. 成都:电子科技大学,2011.

LING Dan. Research on Weibull Distribution and Its Applications in Mechanical Reliability Engineering[D]. Chengdu: University of Electronic Science and Technology of China, 2011.

[7] JOHNSON N L. Systems of Frequency Curves Generated, by Methods of Translation[J]. Biometrika, 1949, 36:149-76.

[8] ZAMAN K, RANGAVAJHALA S, MCDONALD M, et al. A Probabilistic Approach for Representation of Interval Uncertainty[J]. Reliability Engineering and System Safety, 2011, 96(1):117-130.

[9] DEBROTA D J, ROBERTS S D, SWAIN J J, et al. Input Modeling with the Johnson System of Distributions[C]∥Proceedings of the 1988 Winter Simulation Conference. San Diego, 1988:165-179.

[10] 盛骤,谢式千,潘承毅. 概率论与数理统[M].4版. 北京: 高等教育出版社,2008.

SHENG Zhou, XIE Shiqian, PAN Chengyi. Probability Theory & Mathematical Statistics[M]. 4th ed. Beijing: Higher Education Press, 2008.

[11] SLIFKER J F, SHAPIRO S S. The Johnson System: Selection and Parameter Estimation[J]. Technometrics, 1980, 22: 239-246.

[12] 陈开颜,张鹏,邓高明,等. 相关性电磁分析攻击及最小样本量分析[J]. 华中科技大学学报(自然科学版),2011, 39(1):32-35.

CHEN Kaiyan, ZHANG Peng, DENG Gaoming, et al. Electromagnetic Correlation Attacks and Minimum Number and Attack Traces Analyzing[J]. Journal of Huazhong University of Science and Technology (Natural Science Edition), 2011, 39(1):32-35.

[13] 倪延延,张晋昕. 假设检验时样本含量估计中容许误差δ的合理选取[J]. 循证医学,2011, 11(6):370-372.

NI Yanyan, ZHANG Jinxin. How to Determine Permissible Error δ Value Properly When Computing Sample Sizes in Hypothesis Tests[J]. The Journal of Evidence-based Medicine,2011, 11(6):370-372.

[14] 张耀,曹小平. 材料力学[M]. 北京:清华大学出版社,2015.

ZHANG Yao, CAO Xiaoping. Mechanics of Materials[M]. Beijing: Tsinghua University Press, 2015.

Johnson Distribution Fitting Method Based on Minimum Sample Space

ZHANG Qi1 WU Yizhong2 LIU Xin2 QIAO Ping2

1.National NC System Engineering Research Center of Huazhong University of Science & Technology,Wuhan,430074 2.Research Center for National CAD Support Software Engineering and Technology of Huazhong University of Science & Technology,Wuhan,430074

Abstract: When calculating the reliability index based on the stress-strength interference model, statistical methods were used to obtain the probability distribution of stress and strength. However, in the design processes of complex electromechanical products, the acquisition of each sample values was an expensive evaluation, so, the acquisition of a large number of sample values would be confronted with serious problems such as time consuming. Aiming at this problem, the parameterized Johnson distribution family function was introduced to fit the probability distribution, and the Johnson distribution fitting method was studied based on the minimum sample space. The parameters of the Johnson distribution family function were solved by percentile matching method and the minimum sample space was determined based on the hypothesis testing principle. To ensure the fitting accuracy, K-S test was introduced to analyze the consistency between the real distribution and the fitted distribution. The results of engineering applications show that the method may quickly obtain the distribution of samples based on a small number of expensive evaluation samples, which lays a foundation for the rapid calculation of reliability index.

Key words: Johnson distribution family function; hypothesis testing; minimum sample space; conformance analysis

中图分类号TH122

DOI:10.3969/j.issn.1004-132X.2019.21.007

开放科学(资源服务)标识码(OSID):

收稿日期2018-05-15

基金项目国家自然科学基金资助项目(61575205)

(编辑 王艳丽)

作者简介张 琪,男,1994年生,硕士研究生。研究方向为不确定性优化设计。E-mail:hust_zhangqi@hust.edu.cn。吴义忠(通信作者),男,1970年生,教授、博士研究生导师。研究方向为多领域建模与仿真、多学科设计优化、计算机辅助设计与图形学。