基于最小樣本空間的Johnson分布擬合方法

張 琪 吳義忠 劉 鑫 喬 平

1.華中科技大學(xué)國家數(shù)控系統(tǒng)工程技術(shù)研究中心，武漢，4300742.華中科技大學(xué)國家企業(yè)信息化應(yīng)用支撐軟件工程技術(shù)研究中心，武漢，430074

0 引言

現(xiàn)代復(fù)雜機(jī)電產(chǎn)品設(shè)計常受到知識的缺乏和產(chǎn)品運(yùn)行環(huán)境的變化等方面不確定性因素影響，致使設(shè)計的機(jī)電產(chǎn)品在運(yùn)行時部分可靠度指標(biāo)可能會發(fā)生變化或偏移導(dǎo)致引發(fā)產(chǎn)品故障[1]，因此，在產(chǎn)品設(shè)計優(yōu)化階段需基于大量的樣本數(shù)據(jù)來充分考慮各方面的不確定性因素。然而，實際優(yōu)化過程常受到計算資源的制約,在滿足特定條件下應(yīng)盡可能地減少估值采樣的次數(shù)。根據(jù)少量樣本點擬合出總體分布，基于該分布計算可靠度指標(biāo)，可有效避免多次昂貴估值帶來的巨大計算損耗，提高計算效率。

在對總體分布未知的樣本數(shù)據(jù)進(jìn)行分布函數(shù)擬合時，可選擇自適應(yīng)的分布形式對總體進(jìn)行擬合，如標(biāo)準(zhǔn)的二參數(shù)和三參數(shù)Weibull分布，調(diào)整其參數(shù)即可實現(xiàn)分布模型接近于指數(shù)分布、正態(tài)分布等分布模型。在樣本數(shù)據(jù)較多的情況下，圖解法、極大似然估計法[2]、最小二乘法關(guān)于水平殘差和垂直殘差平方和最小[3]等方法對所需參數(shù)的估值較為精確。近年來，Weibull分布在機(jī)電產(chǎn)品的可靠性分析中有著非常重要的應(yīng)用[2-5]，然而Weibull分布中只有兩個或三個參數(shù)，不能全面表達(dá)總體的分布特征，同時，根據(jù)有限的樣本數(shù)據(jù)獲得所需參數(shù)的精確估值還沒有得到有效的解決[6]。

四參數(shù)的Johnson分布族函數(shù)[7](下文簡稱Johnson分布)于1949年被提出，后來發(fā)展成包含四種不同類型的分布族函數(shù)。由于具有多參數(shù)、多類型的特征，Johnson分布通過對樣本擬合可表達(dá)出更多的總體的分布特征。Johnson分布具有較強(qiáng)的自適應(yīng)性，選擇合適的類型和調(diào)整參數(shù)，可接近于任意一個標(biāo)準(zhǔn)連續(xù)型分布模型，其中包括Weibull分布模型[8]。DEBROTA等[9]提出配矩法(moment matching)、百分位數(shù)配比法(percentile matching)、最小二乘法(least squares)和最小范數(shù)估值法(minimum Lpnorm estimation)四種方法，這些方法根據(jù)有限的樣本數(shù)據(jù)即可實現(xiàn)對Johnson分布中四個參數(shù)的精確估計。

本文引用Johnson分布作為未知總體的分布形式對樣本進(jìn)行擬合，并基于假設(shè)檢驗[10]，控制第Ⅰ類錯誤和第Ⅱ類錯誤的發(fā)生概率，根據(jù)Z檢驗法中雙邊檢驗問題的施行特征函數(shù)和曲線，推導(dǎo)出樣本空間的最小值。最后設(shè)計實驗進(jìn)行一致性檢驗，采用K-S(Kolmogorov-Smirnov)擬合檢驗法檢驗總體的擬合分布與真實分布的一致性。實驗結(jié)果表明，當(dāng)樣本空間達(dá)到最小樣本空間時，Johnson分布根據(jù)有限樣本擬合出的總體分布與樣本的真實分布具有一致性，滿足擬合精度要求。

1 Johnson分布族函數(shù)及其擬合

1.1 Johnson分布族函數(shù)

為便于統(tǒng)一表達(dá)各種不同類型的連續(xù)型隨機(jī)變量的累積分布函數(shù)FX(x)=Pr(X≤x)和概率密度函數(shù)f=F′(x)，Johnson提出了四參數(shù)的Johnson分布族函數(shù)。Johnson分布的概率密度函數(shù)為

(1)

其中，γ和δ為形狀參數(shù)；ε為位置參數(shù)；λ為尺度參數(shù)；f(·)為簡單的函數(shù)表達(dá)式，根據(jù)f(·)的不同，Johnson分布可分為對數(shù)正態(tài)(lognormal, SL)、無界 (unbounded, SU)、有界(bounded, SB)和正態(tài)(normal, SN)四種不同類型。

考慮到客觀條件的制約，生產(chǎn)實踐中獲得的數(shù)據(jù)大多數(shù)情況下都服從SB類型的Johnson分布。以SB類型為例，三類參數(shù)對Johnson分布模型的影響如圖1～圖3所示。

圖1 形狀參數(shù)γ、δ對模型的影響Fig.1 The influence of shape parameters γ、δon the model

圖2 位置參數(shù)ε對模型的影響Fig.2 The influence of location parameters εon the model

圖3 尺度參數(shù)λ對模型的影響Fig.3 The influence of dimension parameters λon the model

f(·)和f′(·)的具體表達(dá)式如下：

x的取值范圍H為

對式(1)進(jìn)行積分，即可得到如下各個類型相應(yīng)的累積分布函數(shù)：

無界型累積分布函數(shù)

-∞

正態(tài)型累積分布函數(shù)

-∞

對數(shù)正態(tài)型累積分布函數(shù)

FSL(x)=

有界型累積分布函數(shù)

FSB(x)=

式中，Erf(·)為高斯誤差函數(shù)；Erfc(·)為誤差互補(bǔ)函數(shù)；ArcSinh(·)為反雙曲正弦函數(shù)。

Johnson分布通過變形可得

(2)

該式為Johnson轉(zhuǎn)換式，可將連續(xù)型隨機(jī)變量X映射到一個服從標(biāo)準(zhǔn)正態(tài)分布的隨機(jī)變量Z上[7]。通過Johnson轉(zhuǎn)換式可將樣本數(shù)據(jù){x1,x2,…,xn}轉(zhuǎn)換成一組服從標(biāo)準(zhǔn)正態(tài)分布的新樣本{z1,z2,…,zn}，根據(jù)樣本對總體X進(jìn)行Johnson分布擬合時，擬合精度越高, 則新樣本{z1,z2,…,zn}對標(biāo)準(zhǔn)正態(tài)分布的服從度就越高。

1.2 Johnson分布擬合

對一組樣本的總體分布進(jìn)行擬合時，首先根據(jù)樣本所反映的總體分布特征選擇合適的擬合分布形式，再估算出分布形式中各個參數(shù)的值，從而得到總體的擬合分布。而Johnson分布具有較強(qiáng)的自適應(yīng)性，通過選擇類型和調(diào)整參數(shù)，Johnson分布模型可近似于其他不同形式的連續(xù)型分布模型，即Johnson分布可統(tǒng)一表達(dá)不同形式的連續(xù)型分布。因此，當(dāng)對一組真實分布形式完全未知的樣本進(jìn)行擬合時，不需要分析樣本所反映的總體分布特征來選擇合適的分布形式進(jìn)行總體擬合，直接選用Johnson分布作為總體的分布形式加以擬合即可。

本文通過MATLAB軟件編程實現(xiàn)了Johnson分布對總體擬合的過程。根據(jù)樣本的分布特征，基于Johnson分布類型的選擇準(zhǔn)則正確選擇出合適的分布類型，并得到其相應(yīng)的參數(shù)函數(shù)表達(dá)式，即概率密度函數(shù)表達(dá)式和累積分布函數(shù)表達(dá)式?？紤]到百分位數(shù)配比法比其他方法對所需參數(shù)進(jìn)行估值時更加簡便且能保證估值精度[11]，本文采用百分位數(shù)配比法對所需參數(shù)進(jìn)行精確估計，將各個未知參數(shù)的估值代入?yún)?shù)函數(shù)表達(dá)式中即可得到總體的擬合分布。其擬合流程如圖4所示。

圖4 Johnson分布對總體擬合流程圖Fig.4 The flow chart of Johnson distribution fitting

2 最小樣本空間分析

基于假設(shè)檢驗，通過控制第Ⅰ類錯誤和第Ⅱ類錯誤發(fā)生的概率，根據(jù)Z檢驗法中雙邊檢驗問題的施行特征函數(shù)和曲線，在保證擬合精度的條件下對擬合所需的最小樣本空間進(jìn)行推導(dǎo)。

基于Johnson分布對隨機(jī)變量X的總體分布進(jìn)行擬合后得到的是復(fù)雜的分布函數(shù)表達(dá)式，若直接根據(jù)該分布進(jìn)行假設(shè)檢驗推導(dǎo)出擬合所需最小的樣本空間則更為復(fù)雜。文獻(xiàn)[12]提出一種簡便的解決方法：通過Fisher轉(zhuǎn)換將隨機(jī)變量映射到一個服從正態(tài)分布的隨機(jī)變量上，在正態(tài)分布上推導(dǎo)出擬合所需的最小樣本空間。由于對總體進(jìn)行Johnson分布擬合的精度與通過轉(zhuǎn)換式(2)得到的新樣本對標(biāo)準(zhǔn)正態(tài)分布的服從度相關(guān)，所以本文采用Johnson轉(zhuǎn)換將隨機(jī)變量X映射到一個服從標(biāo)準(zhǔn)正態(tài)分布的隨機(jī)變量Z上，在標(biāo)準(zhǔn)正態(tài)分布上通過控制第Ⅰ類錯誤和第Ⅱ類錯誤發(fā)生的概率得到樣本空間的范圍。為此，引入施行特征函數(shù)[10]。

定義若C是參數(shù)θ的某檢驗問題的一個檢驗法,則稱β(θ)=Pθ(接受H0)為檢驗法C的施行特征函數(shù)或OC函數(shù)，其圖形稱為OC曲線。

根據(jù)上述定義，考慮到雙邊檢驗問題H0：μ=μ0,H1：μ≠μ0, 正態(tài)總體均值的Z檢驗法的OC函數(shù)為

(3)

式中，Φ(·)為求標(biāo)準(zhǔn)正態(tài)分布的累計概率密度函數(shù)。

其OC曲線如圖5所示，β(μ)是|λ|的嚴(yán)格單調(diào)下降函數(shù)。

圖5 OC曲線圖Fig.5 OC curve graph

在雙邊檢驗問題中，若要求對H1中滿足|μ-μ0|≥δ>0的μ處的函數(shù)值β(μ)≤β，則需要解超越方程：

由此知只要樣本空間n滿足

即只要n滿足

(4)

就能使當(dāng)μ∈H1且|μ-μ0|≥δ(δ>0，為取定的值)時，第Ⅰ類錯誤發(fā)生的概率不超過給定的值α，第Ⅱ類錯誤發(fā)生的概率不超過給定的值β。

當(dāng)通過式(4)確定樣本空間范圍時，由文獻(xiàn)[13]可知，容許誤差δ是假設(shè)檢驗所試圖揭示樣本與整體之間的差異大小，總體標(biāo)準(zhǔn)差σ體現(xiàn)個體變異度，對于沒有給定專業(yè)意義上的容許誤差水平的情況，用0.25倍或0.50倍的σ來設(shè)定δ。本文取α=0.05，β=0.05，δ=0.5σ。查表得zα/2=z0.025=1.96，zβ=z0.05=1.645，將上述數(shù)據(jù)代入式(4)得

即通過控制第Ⅰ類錯誤和第Ⅱ類錯誤發(fā)生的概率可得到最小的樣本空間為52。

3 實驗驗證與結(jié)果分析

基于假設(shè)檢驗推導(dǎo)出當(dāng)樣本空間不小于52時，采用Johnson分布對總體的分布進(jìn)行擬合可達(dá)到理論精度?？紤]到樣本空間較小且變量為連續(xù)型時，K-S擬合檢驗法比χ2擬合檢驗法具有更強(qiáng)的檢出力，因此本文采用K-S檢驗法設(shè)計實驗來驗證當(dāng)樣本空間為52時，Johnson分布對總體的擬合是否能夠達(dá)到要求的精度。實驗步驟如下：

(1)采集樣本點。從真實分布已知的總體X中隨機(jī)采集4組樣本, 各組樣本空間分別為42，47，52，57。

(2)確定分布的表達(dá)式。將步驟(1)作為Johnson擬合程序的輸入，運(yùn)行程序, 輸出即為步驟(1)中各組樣本所對應(yīng)的Johnson分布,記擬合分布的總體為X。當(dāng)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檢驗結(jié)果Fig.9 The result of K-S test, if X～N

圖10 X～Γ時，K-S檢驗結(jié)果Fig.10 The result of K-S test, if X～Γ

圖11 X～W時，K-S檢驗結(jié)果Fig.11 The result of K-S test, if X～W

本實驗中，步驟(1)中總體的真實分布形式選取為幾種常見的連續(xù)型分布，包括正態(tài)分布、Gamma分布、Weibull分布等, 實驗結(jié)果如圖9～圖 11所示。根據(jù)圖9～圖11所示的實驗檢驗結(jié)果，通過分析可知:

(1)不失一般性，隨著樣本空間增大，擬合精度提高，K-S檢驗法的檢驗結(jié)果q值呈增大趨勢。當(dāng)樣本空間為52時，檢驗結(jié)果q值均大于0.05，說明假設(shè)成立，從X和X中分別隨機(jī)抽取的兩組樣本來自于同一分布。這表明，樣本空間為52時，總體的擬合分布與真實分布具有一致性，即在一定精度條件下，擬合分布等價于真實分布。

(2)當(dāng)樣本空間為52，樣本數(shù)據(jù)取自真實分布X～N(2,2)和X～W(3,5)時，K-S檢驗法得到的q值遠(yuǎn)大于0.05，說明實際擬合精度遠(yuǎn)高于要求精度。但一般情況下，q值均以0.05為下界，在一定范圍內(nèi)波動。這表明，采用Johnson分布對一組總體分布形式完全未知的樣本進(jìn)行擬合時，52作為最小樣本空間具有普適性。

4 工程算例

曲軸是發(fā)動機(jī)中最重要的部件，它的可靠性對發(fā)動機(jī)的使用壽命有著重要的影響。當(dāng)基于應(yīng)力-強(qiáng)度干涉模型計算曲柄頸的可靠度指標(biāo)時，需在工作過程中統(tǒng)計曲柄頸危險截面處應(yīng)力值的變化從而獲得其概率密度。將曲柄臂抽象為矩形，如圖12所示，當(dāng)曲軸各參數(shù)取表1中的值時，已知曲柄頸處所受的力F服從正態(tài)分布N(15,1)，即可根據(jù)上述數(shù)據(jù)構(gòu)建曲軸的力學(xué)模型，進(jìn)而求得所需應(yīng)力值。

圖12 曲軸力學(xué)模型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

利用測力傳感器，從服從正態(tài)分布N(15,1)的F中隨機(jī)抽取52個值作為樣本數(shù)據(jù)點，通過力學(xué)分析和公式計算[14]，曲柄頸危險截面處對應(yīng)的52個應(yīng)力值如表2所示。

采用Johnson分布對所求的52個曲柄頸處所受力F進(jìn)行擬合，得到應(yīng)力F的概率密度函數(shù)

采用Johnson分布對所求的52個應(yīng)力值進(jìn)行擬合，得到應(yīng)力σ0的概率密度函數(shù)為

表2 52個F樣本點及其對應(yīng)的危險截面處應(yīng)力σ0

力F和應(yīng)力σ0的概率分布擬合如圖13和圖14所示。

圖13 力F的概率分布擬合Fig.13 The fitting pdf of F

圖14 應(yīng)力σ0的概率分布擬合Fig.14 The fitting pdf of σ0

5 結(jié)論

通過引入Johnson分布函數(shù)，采用百分位數(shù)配比法求解Johnson分布函數(shù)的參數(shù)，并基于假設(shè)檢驗原理確定滿足精度要求的最小樣本空間，同時通過K-S擬合檢驗法對總體真實分布與擬合分布進(jìn)行一致性驗證，保證擬合精度滿足要求，實現(xiàn)了Johnson分布對總體分布擬合過程。本文的研究工作可快速確定對總體分布擬合所需的最小樣本空間，提高基于少量昂貴估值對總體分布擬合的精度和效率，為可靠性設(shè)計優(yōu)化過程中的可靠度計算提供理論依據(jù)，對提高可靠性設(shè)計優(yōu)化的計算效率具有重要意義。