王 丹，李衛(wèi)國

(北京航空航天大學(xué)數(shù)學(xué)與系統(tǒng)科學(xué)學(xué)院，北京 100191)

在航空領(lǐng)域中，如對導(dǎo)彈推進(jìn)劑對射程的影響和飛行器失效率的研究，在火工品領(lǐng)域中，如對軍火用品失效率隨時(shí)間變化的研究，由于實(shí)際環(huán)境與成本等原因，不適宜進(jìn)行大量的試驗(yàn)，試驗(yàn)次數(shù)受到限制，繼而生成小樣本。在進(jìn)行必要的統(tǒng)計(jì)推斷過程中，參數(shù)估計(jì)是很重要的一部分，通常我們需要用一組已知的樣本 x1，x2，L，xn，來估計(jì)總體X的參數(shù)θ，估計(jì)方法有很多種，比如最大似然估計(jì)、最小二乘估計(jì)或者矩估計(jì)等。如果試驗(yàn)的次數(shù)足夠多，得到的樣本足夠大，那么根據(jù)大數(shù)定律，樣本的估計(jì)值將收斂于真正的參數(shù)值，但是，由于試驗(yàn)次數(shù)的限制，得到的用于進(jìn)行估計(jì)參數(shù)的樣本就比較有限，那么根據(jù)隨機(jī)性，推斷得到的結(jié)果很可能與實(shí)際情況完全不符合。比如，在正常情況下，飛行器失效率應(yīng)該隨時(shí)間遞增，但根據(jù)試驗(yàn)數(shù)據(jù)所得到的失效率卻并不滿足遞增的規(guī)律，在這種情況下，應(yīng)該盡可能提取樣本所含信息量，并將其用于統(tǒng)計(jì)推斷。此時(shí)，我們自然會(huì)考慮是否可以調(diào)整數(shù)據(jù)使之滿足實(shí)際中應(yīng)該符合的關(guān)系，而保序回歸算法——PAVA(Pool Adjacent Violators Algorithm)算法可以用來調(diào)整數(shù)據(jù)。

保序回歸的研究是從20世紀(jì)中葉在國外開始的，1955年Ayer等人提出了關(guān)于保序回歸的求解算法——PAVA算法，后來保序回歸問題得到廣泛的重視。1972年由Brunk，H.D，Bartholomew，D.J等人聯(lián)合編著了《The theory and application of Isotonic Regression》［1］，這本書總結(jié)了之前所有關(guān)于保序問題的研究，得到很熱烈的反響。1988年由Tim Robertson，F(xiàn).T.Wright和R.L.Dykstra聯(lián)合編著的《Oder Restricted Statistical Inference》使保序回歸有了進(jìn)一步發(fā)展。在國內(nèi)，東北師范大學(xué)的史寧中教授在1993年02期《應(yīng)用概率統(tǒng)計(jì)》中發(fā)表的文章《保序回歸與最大似然估計(jì)》［2］對保序回歸做了比較完整的介紹，他用實(shí)例導(dǎo)出統(tǒng)計(jì)模型，系統(tǒng)地總結(jié)了保序回歸的性質(zhì)、求解方法和其與最大似然估計(jì)之間的關(guān)系，并把問題擴(kuò)展到多維保序回歸與廣義保序回歸?，F(xiàn)在仍有許多人致力于保序問題的研究，并且可以看到，保序回歸的思想已經(jīng)應(yīng)用到各個(gè)科學(xué)領(lǐng)域。

1 保序回歸PAVA算法簡介

變量滿足序關(guān)系，但觀察或經(jīng)過統(tǒng)計(jì)推斷得到的數(shù)據(jù)不滿足本該有的順序時(shí)，對數(shù)據(jù)通過加權(quán)平均的方法進(jìn)行調(diào)整，這就是保序回歸PAVA算法所解決的問題。例如u1，u2，…，uk是一組變量，且滿足u1≤u2≤…≤uk，而是u1，u2，L，uk估計(jì)值，那么利用 PAVA 算法進(jìn)行保序回歸的過程如下:

(3)由此可以得到變量 u1，u1，…，uk經(jīng)PAVA算法調(diào)整后的值。

保序回歸解的唯一性已被驗(yàn)證［3］，上述過程得到的u1，u2，…，uk滿足順序關(guān)系的解是唯一的。

下面給出一個(gè)有關(guān)飛行器失效率的例子，為測定某飛行器失效率隨時(shí)間變化的規(guī)律，一般是在給定的時(shí)間水平 t1，t2，…，tk上，并且 t1＜t2＜…＜ tk，測試［0，ti］時(shí)間產(chǎn)品失效的率，所得數(shù)據(jù)如下:

表1 飛行器產(chǎn)品失效系列試驗(yàn)數(shù)據(jù)表

產(chǎn)品的平均失效率一般會(huì)隨著時(shí)間段的長度的增加而遞增，但根據(jù)此次是試驗(yàn)結(jié)果，估計(jì)的失效率并沒有呈現(xiàn)出單調(diào)遞增的趨勢，此時(shí)我們采用PAVA算法進(jìn)行調(diào)整。

在應(yīng)用PAVA算法時(shí)，如果數(shù)據(jù)量比較大，應(yīng)該利用編程解決，由于樣本的不確定，根據(jù)定義對PAVA算法編程就變得非常困難，此時(shí)就可以用有約束線性最小二乘函數(shù)lsqlin來實(shí)現(xiàn)［3］。

2 問題的提出

對數(shù)據(jù)進(jìn)行調(diào)整之后再作分析的過程中，如對火工品進(jìn)行分析時(shí)，假設(shè)火工品在其貯存時(shí)間t時(shí)刻產(chǎn)品的臨界刺激量為x(t)，且x(t)服從正態(tài)分布 N(μ(t)，σ2(t))，在工程應(yīng)用時(shí)，假設(shè) μ(t)和σ2(t)分別關(guān)于f(t)和g(t)是線性模型，其中f(t)和g(t)是關(guān)于t的已知函數(shù)，一般可由產(chǎn)品的歷史信息或者相似產(chǎn)品信息確定，接下來再考慮可靠性，模型如下［4］:

此處就出現(xiàn)滿足序關(guān)系的變量關(guān)于自變量的線性回歸，有一個(gè)問題自然需要考慮，就是在線性回歸中，因變量是估計(jì)值直接關(guān)于自變量回歸呢，還是經(jīng)過PAVA算法調(diào)整之后再回歸，兩種結(jié)果一致嗎?

更一般地，由線性回歸的定義可以知道線性回歸的因變量隨著自變量的增加一定是滿足遞增或遞減的關(guān)系，那么在這種情況下應(yīng)用最小二乘線性回歸時(shí)，先對觀察到的因變量使用PAVA算法進(jìn)行調(diào)整是否有意義呢?

接下來就應(yīng)用數(shù)學(xué)實(shí)驗(yàn)進(jìn)行模擬，以期望得出結(jié)論。

3 實(shí)驗(yàn)?zāi)M及結(jié)論

在模擬之前，首先假定因變量Y與因變量X滿足關(guān)系:Y=X+1，試驗(yàn)樣本為100，具體步驟如下:

(1)為了要生成100個(gè)隨機(jī)樣本，首先使用軟件生成100個(gè)服從正態(tài)分布的隨機(jī)數(shù);

(2)取100 個(gè)自變量 0.1，0.2，0.4，…，10，相應(yīng)能得到 100 個(gè)因變量 1.1，1.2，1.3，…，11，現(xiàn)在將(1)中產(chǎn)生的100個(gè)服從正態(tài)分布的隨機(jī)數(shù)加在得到因變量中，就能得到第一組的樣本?？紤]到自變量間隔不同，所加的隨機(jī)部分對因變量順序影響也不同，為了說明問題，對自變量間隔從0.1，0.2一直到0.9這9中情況都作考慮;

(3)對于(2)中得到的每組數(shù)據(jù)，直接進(jìn)行最小二乘回歸，然后用得到的模型所對應(yīng)的值與真實(shí)值進(jìn)行比較，計(jì)算出偏差絕對值之和;

(4)對于(3)中得到的每組數(shù)據(jù)，首先使用PAVA算法進(jìn)行調(diào)整，然后用調(diào)整過的值進(jìn)行最小二乘回歸，再比較模型預(yù)測值與真實(shí)值，計(jì)算出偏差絕對值;

(5)重復(fù)步驟(1)至(4)多次，對實(shí)驗(yàn)中絕對值之和分別進(jìn)行比較，得出結(jié)論。

下面展示其中一組試驗(yàn)的數(shù)據(jù):

在這組數(shù)據(jù)中，random列為隨機(jī)抽取的服從正態(tài)分布的數(shù)列，y10為Y的真實(shí)值，y11相當(dāng)于隨機(jī)抽取的樣本，y12是將y11經(jīng)過PAVA算法調(diào)整之后得到的數(shù)據(jù)，PER_11是y11關(guān)于x1回歸所得模型的預(yù)測值，PRE_12是y12關(guān)于x1回歸所得模型的預(yù)測值，shengyu11和shengyr12分別為預(yù)測值與真實(shí)值之間差的絕對值，可以看到兩中不同的處理方法中，偏差絕對值之和分別被4.69和4.689 41。

表2 部分試驗(yàn)數(shù)據(jù)表

4 結(jié)論

為保證實(shí)驗(yàn)結(jié)果的準(zhǔn)確性，我們可以多做幾次模擬，比較采用兩種不同處理方法時(shí)，預(yù)測值與真實(shí)值偏差絕對值之和，以期望得到結(jié)論。在本實(shí)驗(yàn)中，我們做了10次模擬，并將部分結(jié)果展示如下:

由上述結(jié)果我們得出如下結(jié)論，兩種方法的擬合效果是不分上下的，在實(shí)際應(yīng)用中，如果數(shù)據(jù)量比較大，使用PAVA算法作調(diào)整并不能將擬合效果提升反而會(huì)增加相當(dāng)多的工作量，所以建議在因變量滿足序關(guān)系的情況下作最小二乘回歸時(shí)，無需對因變量使用PAVA算法作調(diào)整，應(yīng)直接做最小二乘回歸。

表3 試驗(yàn)結(jié)果對比表

［1］ Barlow R E，Bartholotemew D J，Bremner J M，et al.Statitical Inference under Order Restrictions:The Theory and Application of Isotonic Regression［M］.New York:Wilely，1972.

［2］史寧中.保序回歸與最大似然估計(jì)［J］.應(yīng)用概率統(tǒng)計(jì)，1993(2):203－215.

［3］邢務(wù)強(qiáng).保序回歸的研究及應(yīng)用［D］.西安:西北工業(yè)大學(xué)，2002.

［4］趙鳳治，尉繼英.約束最優(yōu)化計(jì)算方法［M］.北京:北京科學(xué)出版社，1991.

［5］ Shi Ningzhong.Maximum likelihood estimation of means and variances from normal populations under simultaneous order restrictions［J］.Journal of Multivariate Analysis，1996，50(2):282 －293.

［6］ Shi Ningzhong，Jiang Hua.Maximum likelihood estimation of isotonic normal means with unknown variances［J］.Journal of Multivariate Analysis，1998，64(2):183－193.

［7］張潤楚.多元統(tǒng)計(jì)分析［M］.北京:北京科學(xué)出版社，2006.

［8］洪東跑，趙宇，溫玉全.基于序約束的火工品可靠性試驗(yàn)數(shù)據(jù)分析［J］.含能材料，2008，16(5):556－559.

［9］洪東跑，趙宇，溫玉全.利用感度試驗(yàn)數(shù)據(jù)分析火工品貯存可靠性［J］.含能材料，2009，17(5):608－611.