基于最小樣本空間的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分布對總體分布擬合過程。本文的研究工作可快速確定對總體分布擬合所需的最小樣本空間，提高基于少量昂貴估值對總體分布擬合的精度和效率，為可靠性設計優化過程中的可靠度計算提供理論依據，對提高可靠性設計優化的計算效率具有重要意義。