保序调整对线性回归影响的试验分析

王 丹，李卫国

(北京航空航天大学数学与系统科学学院，北京 100191)

在航空领域中，如对导弹推进剂对射程的影响和飞行器失效率的研究，在火工品领域中，如对军火用品失效率随时间变化的研究，由于实际环境与成本等原因，不适宜进行大量的试验，试验次数受到限制，继而生成小样本。在进行必要的统计推断过程中，参数估计是很重要的一部分，通常我们需要用一组已知的样本 x1，x2，L，xn，来估计总体X的参数θ，估计方法有很多种，比如最大似然估计、最小二乘估计或者矩估计等。如果试验的次数足够多，得到的样本足够大，那么根据大数定律，样本的估计值将收敛于真正的参数值，但是，由于试验次数的限制，得到的用于进行估计参数的样本就比较有限，那么根据随机性，推断得到的结果很可能与实际情况完全不符合。比如，在正常情况下，飞行器失效率应该随时间递增，但根据试验数据所得到的失效率却并不满足递增的规律，在这种情况下，应该尽可能提取样本所含信息量，并将其用于统计推断。此时，我们自然会考虑是否可以调整数据使之满足实际中应该符合的关系，而保序回归算法——PAVA(Pool Adjacent Violators Algorithm)算法可以用来调整数据。

保序回归的研究是从20世纪中叶在国外开始的，1955年Ayer等人提出了关于保序回归的求解算法——PAVA算法，后来保序回归问题得到广泛的重视。1972年由Brunk，H.D，Bartholomew，D.J等人联合编著了《The theory and application of Isotonic Regression》［1］，这本书总结了之前所有关于保序问题的研究，得到很热烈的反响。1988年由Tim Robertson，F.T.Wright和R.L.Dykstra联合编著的《Oder Restricted Statistical Inference》使保序回归有了进一步发展。在国内，东北师范大学的史宁中教授在1993年02期《应用概率统计》中发表的文章《保序回归与最大似然估计》［2］对保序回归做了比较完整的介绍，他用实例导出统计模型，系统地总结了保序回归的性质、求解方法和其与最大似然估计之间的关系，并把问题扩展到多维保序回归与广义保序回归。现在仍有许多人致力于保序问题的研究，并且可以看到，保序回归的思想已经应用到各个科学领域。

1 保序回归PAVA算法简介

变量满足序关系，但观察或经过统计推断得到的数据不满足本该有的顺序时，对数据通过加权平均的方法进行调整，这就是保序回归PAVA算法所解决的问题。例如u1，u2，…，uk是一组变量，且满足u1≤u2≤…≤uk，而是u1，u2，L，uk估计值，那么利用 PAVA 算法进行保序回归的过程如下:

(3)由此可以得到变量 u1，u1，…，uk经PAVA算法调整后的值。

保序回归解的唯一性已被验证［3］，上述过程得到的u1，u2，…，uk满足顺序关系的解是唯一的。

下面给出一个有关飞行器失效率的例子，为测定某飞行器失效率随时间变化的规律，一般是在给定的时间水平 t1，t2，…，tk上，并且 t1＜t2＜…＜ tk，测试［0，ti］时间产品失效的率，所得数据如下:

表1 飞行器产品失效系列试验数据表

产品的平均失效率一般会随着时间段的长度的增加而递增，但根据此次是试验结果，估计的失效率并没有呈现出单调递增的趋势，此时我们采用PAVA算法进行调整。

在应用PAVA算法时，如果数据量比较大，应该利用编程解决，由于样本的不确定，根据定义对PAVA算法编程就变得非常困难，此时就可以用有约束线性最小二乘函数lsqlin来实现［3］。

2 问题的提出

对数据进行调整之后再作分析的过程中，如对火工品进行分析时，假设火工品在其贮存时间t时刻产品的临界刺激量为x(t)，且x(t)服从正态分布 N(μ(t)，σ2(t))，在工程应用时，假设 μ(t)和σ2(t)分别关于f(t)和g(t)是线性模型，其中f(t)和g(t)是关于t的已知函数，一般可由产品的历史信息或者相似产品信息确定，接下来再考虑可靠性，模型如下［4］:

此处就出现满足序关系的变量关于自变量的线性回归，有一个问题自然需要考虑，就是在线性回归中，因变量是估计值直接关于自变量回归呢，还是经过PAVA算法调整之后再回归，两种结果一致吗?

更一般地，由线性回归的定义可以知道线性回归的因变量随着自变量的增加一定是满足递增或递减的关系，那么在这种情况下应用最小二乘线性回归时，先对观察到的因变量使用PAVA算法进行调整是否有意义呢?

接下来就应用数学实验进行模拟，以期望得出结论。

3 实验模拟及结论

在模拟之前，首先假定因变量Y与因变量X满足关系:Y=X+1，试验样本为100，具体步骤如下:

(1)为了要生成100个随机样本，首先使用软件生成100个服从正态分布的随机数;

(2)取100 个自变量 0.1，0.2，0.4，…，10，相应能得到 100 个因变量 1.1，1.2，1.3，…，11，现在将(1)中产生的100个服从正态分布的随机数加在得到因变量中，就能得到第一组的样本。考虑到自变量间隔不同，所加的随机部分对因变量顺序影响也不同，为了说明问题，对自变量间隔从0.1，0.2一直到0.9这9中情况都作考虑;

(3)对于(2)中得到的每组数据，直接进行最小二乘回归，然后用得到的模型所对应的值与真实值进行比较，计算出偏差绝对值之和;

(4)对于(3)中得到的每组数据，首先使用PAVA算法进行调整，然后用调整过的值进行最小二乘回归，再比较模型预测值与真实值，计算出偏差绝对值;

(5)重复步骤(1)至(4)多次，对实验中绝对值之和分别进行比较，得出结论。

下面展示其中一组试验的数据:

在这组数据中，random列为随机抽取的服从正态分布的数列，y10为Y的真实值，y11相当于随机抽取的样本，y12是将y11经过PAVA算法调整之后得到的数据，PER_11是y11关于x1回归所得模型的预测值，PRE_12是y12关于x1回归所得模型的预测值，shengyu11和shengyr12分别为预测值与真实值之间差的绝对值，可以看到两中不同的处理方法中，偏差绝对值之和分别被4.69和4.689 41。

表2 部分试验数据表

4 结论

为保证实验结果的准确性，我们可以多做几次模拟，比较采用两种不同处理方法时，预测值与真实值偏差绝对值之和，以期望得到结论。在本实验中，我们做了10次模拟，并将部分结果展示如下:

由上述结果我们得出如下结论，两种方法的拟合效果是不分上下的，在实际应用中，如果数据量比较大，使用PAVA算法作调整并不能将拟合效果提升反而会增加相当多的工作量，所以建议在因变量满足序关系的情况下作最小二乘回归时，无需对因变量使用PAVA算法作调整，应直接做最小二乘回归。

表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］.应用概率统计，1993(2):203－215.

［3］邢务强.保序回归的研究及应用［D］.西安:西北工业大学，2002.

［4］赵凤治，尉继英.约束最优化计算方法［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］张润楚.多元统计分析［M］.北京:北京科学出版社，2006.

［8］洪东跑，赵宇，温玉全.基于序约束的火工品可靠性试验数据分析［J］.含能材料，2008，16(5):556－559.

［9］洪东跑，赵宇，温玉全.利用感度试验数据分析火工品贮存可靠性［J］.含能材料，2009，17(5):608－611.