基于循序I型刪失數(shù)據(jù)的廣義Pareto分布最優(yōu)刪失計(jì)劃

程從華,程麗娟

(1. 肇慶學(xué)院 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，廣東 肇慶 526061；2. 嶺南師范學(xué)院 數(shù)學(xué)與統(tǒng)計(jì)學(xué)院，廣東 湛江524048)

廣義Paretto分布最早由Pickands提出[1].隨機(jī)變量X服從廣義Pareto分布,如果它的概率密度函數(shù)(PDF)為

(1)

f(x;α,λ)=αλ(1+λx)-(α+1).

(2)

對(duì)應(yīng)的分布函數(shù)為：

F(x;α,λ)=1-(1+λx)-α,x,α,λ>0.

(3)

廣義Pareto模型在極端事件分析中有著廣泛的應(yīng)用.比如，保險(xiǎn)分析中的大額報(bào)單索賠問(wèn)題以及可靠性分析中的失效時(shí)間建模問(wèn)題都可以利用廣義Pareto模型來(lái)進(jìn)行建模分析. Harris研究了保修服務(wù)時(shí)間決策問(wèn)題[2]. Davis和Feldstein利用廣義Pareto模型研究了循序刪失情形下的等效元件失效時(shí)間問(wèn)題[3]. Hosking和Wallis研究了廣義Pareto模型的參數(shù)和分位數(shù)估計(jì)問(wèn)題[4]. Smith研究了非正則條件下的分布族參數(shù)最大似然估計(jì)問(wèn)題[5]. Liang基于非參數(shù)經(jīng)驗(yàn)貝葉斯方法討論了廣義Pareto分布的尺度參數(shù)估計(jì)問(wèn)題[6]. Nigm等在兩樣本和隨機(jī)樣本容量條件下討論了廣義Pareto分布未知參數(shù)的貝葉斯區(qū)間估計(jì)問(wèn)題[7].Wu等在循序刪失數(shù)據(jù)情形下研究了廣義Pareto分布參數(shù)的區(qū)間估計(jì)問(wèn)題[8].

尋找最優(yōu)刪失計(jì)劃是一個(gè)近年來(lái)受到廣泛關(guān)注的問(wèn)題.在循序刪失情形下，本文首先討論廣義Pareto分布未知參數(shù)的最大似然估計(jì)問(wèn)題，并且給出廣義Pareto模型在循序刪失條件下的期望Fisher信息矩陣.其次利用期望Fisher信息矩陣，在三種不同準(zhǔn)則下，討論最優(yōu)刪失計(jì)劃的設(shè)計(jì)問(wèn)題.在壽命分析實(shí)驗(yàn)中，大多數(shù)壽命實(shí)驗(yàn)還會(huì)受到實(shí)驗(yàn)經(jīng)費(fèi)預(yù)算的約束.近幾年來(lái)許多學(xué)者也對(duì)此進(jìn)行了研究，比如：Tse等[9]，Chen等[10]和 Wu等[11]. 在實(shí)驗(yàn)經(jīng)費(fèi)不超過(guò)給定數(shù)額條件下，討論循序刪失計(jì)劃的最優(yōu)設(shè)計(jì)問(wèn)題. 針對(duì)循序I型區(qū)間刪失特點(diǎn)，主要考慮三個(gè)問(wèn)題,分別是如何確定參與壽命分析實(shí)驗(yàn)的元件個(gè)數(shù)，觀測(cè)區(qū)間個(gè)數(shù)以及實(shí)驗(yàn)檢測(cè)區(qū)間長(zhǎng)度.最后給出完成壽命測(cè)試實(shí)驗(yàn)的一個(gè)具體算法，并且通過(guò)一個(gè)具體實(shí)例來(lái)演示本文的方法.

1 期望信息矩陣

為了計(jì)算未知參數(shù)的Fisher信息，需要以下的一些預(yù)備知識(shí)[12].

Xi|Xi-1,Xi-2,…，X1，Ri-1,Ri-2,…，R1～B(Mi,qi),

E(M1)=n,E(R1)=np1(1-q1),E(Xi)=E(Mi)qi,i=1,2,…，m,

利用以上結(jié)論，可以得到期望Fisher信息矩陣E，且E可以表示為：

(4)

這里

(5)

2 最優(yōu)刪失計(jì)劃設(shè)計(jì)

2.1 算法設(shè)計(jì)

在實(shí)驗(yàn)成本約束條件下，這一小節(jié)討論最優(yōu)的實(shí)驗(yàn)設(shè)計(jì)問(wèn)題.利用Ng等定義的如下三個(gè)最優(yōu)準(zhǔn)則來(lái)進(jìn)行實(shí)驗(yàn)設(shè)計(jì)[13].

(2)T-最優(yōu)：最小化協(xié)方差矩陣的跡，tr(V(α,λ))=V11+V22.

(3)F-最優(yōu)：最大化參數(shù)最大似然估計(jì)期望Fisher矩陣的跡，tr(E(α,λ))=E11+E22.

假定相鄰的觀測(cè)區(qū)間長(zhǎng)度差都等于給定的長(zhǎng)度t，同時(shí)假定有n個(gè)元件投入壽命測(cè)試實(shí)驗(yàn)，有m個(gè)觀測(cè)時(shí)刻點(diǎn)，第i個(gè)觀測(cè)區(qū)間的時(shí)間長(zhǎng)度為it，i=1,2,…，m.同時(shí)假設(shè)以下實(shí)驗(yàn)設(shè)計(jì)參數(shù).

(a) 樣本成本：令Cs是每一個(gè)參與測(cè)試的元件價(jià)格，則樣本總成本為nCs.

(b) 檢測(cè)成本：令Ci是每一個(gè)參與測(cè)試的元件檢測(cè)成本，則總檢測(cè)成本為mCi.

因此，壽命測(cè)試實(shí)驗(yàn)的總成本是：

顯而易見(jiàn)的是，每個(gè)實(shí)驗(yàn)準(zhǔn)則都是n,m,t的函數(shù)，記為G(n,m,t).當(dāng)實(shí)驗(yàn)總成本是給定參數(shù)Cb時(shí)，則約束條件變?yōu)椋?/p>

(6)

因此，這個(gè)實(shí)驗(yàn)的最優(yōu)設(shè)計(jì)可以表述為：

(7)

其中，N是正整數(shù).

可以看到目標(biāo)函數(shù)和約束條件都是非線性函數(shù). 下面將利用非線性混合規(guī)劃方法求解上述目標(biāo)函數(shù).非線性規(guī)劃問(wèn)題由Kamat和 Mesquita首先提出[14].關(guān)于非線性混合規(guī)劃方法比較全面的知識(shí)可以參考Grossmann的介紹[15]. 對(duì)本文涉及的目標(biāo)函數(shù)和具體問(wèn)題,主要參考Taha的方法[16]. 基于上述介紹，給出如下算法.

(Ⅳ)對(duì)于給定的n，計(jì)算函數(shù)G(n,m,tmn)的值.

(Ⅶ)計(jì)算最優(yōu)函數(shù)G(n,m,tmn)的值，這里的

則(n*,m*,t*)就是我們尋找的最優(yōu)刪失計(jì)劃.

2.2 演示實(shí)例

min imize(max imize)G(n,m,t),

subject to 10n+5m+m(m+1)t≤600,n,m∈N,andt>0.

利用上節(jié)給出的算法，可以獲得最優(yōu)刪失計(jì)劃如下.

D-最優(yōu)：n*=55,m*=9,t*=0.0556.

T-最優(yōu)：n*=54,m*=11,t*=0.0379.

F-最優(yōu)：n*=56,m*=6,t*=0.2381.

通過(guò)數(shù)值實(shí)驗(yàn)結(jié)果，可以發(fā)現(xiàn)無(wú)論是D-最優(yōu)，T-最優(yōu)還是F-最優(yōu)，本文給出的方法均是可以實(shí)現(xiàn)的.但各方案實(shí)現(xiàn)的具體結(jié)果有所差異，其中D-最優(yōu)方案和T-最優(yōu)方案在結(jié)果上更為相近，F(xiàn)-最優(yōu)方案則有較大差異.這一結(jié)果并不令人意外，因?yàn)镕-最優(yōu)方案利用的是期望Fisher信息矩陣，而D-最優(yōu)方案和T-最優(yōu)方案使用同一個(gè)協(xié)方差矩陣.因此在實(shí)踐中，在小樣本情形時(shí)，建議使用D-最優(yōu)方案和T-最優(yōu)方案,反之使用F-最優(yōu)方案.

[1] Pickands J. Statistical inference using extreme order statistics[J]. The Annals of Statistics, 1975, 3(1): 119-131.

[2] Harris C M. The Pareto distribution as a queue discipline[J]. Operations Research, 1968.16(2): 307-313.

[3] Davis H T, Feldstein M L. The generalized Pareto law as a model for progroressively censored survival data[J]. Biometrika, 1979, 66(2): 299-306.

[4] Hosking J R M,Wallis J R. Parameter and quantile estimation for the generalized Pareto distribution[J]. Technometrics, 1987, 29(3): 339-349.

[5] Smith R L. Maximum likelihood estimation in a class of nonregular cases[J].Bimetrika, 1985, 72(1): 67-90.

[6] Liang T C.Convergence rates for empirical Bayes estimation of the scaleparameter in a Pareto distribution[J]. Computational Statistics & Data Analysis , 1993, 16(1): 35-45.

[7] Nigm A M, Al-Hussaini E K, Jaheen Z F. Bayesian two-sample predictionunder the Lomax model with fixed and random sample size[J]. Journal of Applied Statistics, 2003, 37(6): 527-536.

[8] Wu S F. Interval estimation for the Pareto distribution based on the progressive type II censored sample[J]. Journal of Statistical Computation and Simulation, 2010, 80(4):463-474.

[9] Tse S K, Yang C, Yuen H K. Design and analysis of survival data underan integrated type II interval censoring scheme[J]. Journal of Biopharmaceutical Statistics, 2002, 12(3): 333-345.

[10] Chen J W, Li K H, Lam Y. Bayesian single and double variable sampling plans forthe Weibull distribution with censoring[J]. European Journal of Operational Research, 2007, 177(2): 1062-1073.

[11] Wu S J, Huang S R. Optimal progressive group-censoring pans for exponential distribution in presence of cost constraint[J]. Statistical Papers, 2010, 51(2):4 31-443.

[12] Aggarwala R. Progressive interval censoring: some mathematical results with applications to inference[J]. Communications in Statistics-Theory and Methods, 2001, 30(8): 1921-1935.

[13] Ng H K T, Chan P S, Balakrishnan N. Optimal progressive censoring plan for the Weibull distribution[J]. Technometrics, 2004, 46(6): 470-481.

[14] Kamat M P, Mesquita L.Nonlinear mixed integer programming[M]// Adeli H.Advanced in design optimization. London：Chapman & Hall Press, 1994.

[15] Grossmann I E. Review of nonlinear mixed-integer and disjunctive programming techniques[J]. Optimization and Engineering, 1965, 3(3): 227-252.

[16] Taha H A. Operations research: an introduction[M]. 5th ed.New York:Macmillan Press, 1992.

[17] Dempster A P, Laird N M, Rubin D B. Maximum likelihood from incomplete data via the EM algorithm[J]. Journal of the Royal Statistical Society：Series B, 1977, 39(1): 1- 38.