王 丹，李衛(wèi)國

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

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

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

1 保序回歸PAVA算法簡介

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

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

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

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

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

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

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

2 問題的提出

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

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

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

接下來就應用數(shù)學實驗進行模擬，以期望得出結論。

3 實驗模擬及結論

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

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

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

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

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

(5)重復步驟(1)至(4)多次，對實驗中絕對值之和分別進行比較，得出結論。

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

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

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

4 結論

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

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

表3 試驗結果對比表

［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］.應用概率統(tǒng)計，1993(2):203－215.

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

［4］趙鳳治，尉繼英.約束最優(yōu)化計算方法［M］.北京:北京科學出版社，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)計分析［M］.北京:北京科學出版社，2006.

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

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