定时截尾数据回归分析方法

文章编号:100028055 (2010) 0120142206
定时截尾数据回归分析方法
傅惠民, 岳晓蕊
(北京航空航天大学小样本技术研究中心, 北京100191)
摘要: 提出定时截尾数据回归分析方法,建立定时截尾数据回归方程,给出回归系数和标准差的最佳无偏整体估计、百分位值的点估计及其置信限估计. 详细讨论了工程中常见的极值分布、Wei b ull 分布、正态分布以及一般位置2尺度分布场合下定时截尾数据的回归分析问题. 该方法可以充分利用截尾时刻的试验信息, 并将各个状态下的定时截尾数据作为一个整体进行统计分析,具有信息量大、精度高的特点,从而可以在试验时间、经费等有限的情况下对产品进行高精度的寿命预测和可靠性评估.
关键词: 定时截尾数据; 回归分析; 极值分布; Wei bull 分布; 正态分布; 位置2尺度分
布;
可靠性评估
中图分类号: TB301 ; O211 文献标识码: A
R egression anal y sis method f or type2 Ⅰcensored d ata
FU H ui2mi n , YU E Xiao2r ui
( Re s ea r ch Ce n t e r of Small Sa mp l e Tec h n olo g y ,
Beiji n g U n ive r s it y of A e r o n a u tic s a n d A s t r o n a u tic s , Beiji n g 100191 , Chi n a) Abstract : Re g re s sio n a n al y si s met h o d fo r t y p e2 Ⅰce n s o r ed dat a wa s p r e s e n t e d. The re2 gre s sio n equatio n fo r t yp e2 Ⅰce n so red dat a wa s e s t a bli s hed. The be st unbia sed i nt e g ral e s ti2 mato r s of t he regre s sio n co e fficie nt s a n d t he st a nda r d deviatio n , t he poi nt e sti mato r s a n d co nfi de nce i nt e r val s of p erce ntile s we re give n . The regre s sio n a nal y s i s met h o d s fo r t y p e2Ⅰ ce n so red dat a f ro m e xt re me di st ri butio n , Wei bull di st ri butio n , no r mal di st ri b utio n a n d l oca2 tio n2scale di s t r i b utio n we r e di s c u s sed i n det a il . A s t h e e xp e r i me n t a l i n fo r matio n at t h e ce n s o2 ri ng ti me a nd t he cro s s i nfo r matio n a mo ng t e s t dat a i n diff ere nt st at e s we re f u ll y e x p l oit e d , t h e p re se nt e d met h o d co ul d reach a hi gh p reci sio n . The met h o d p ro v i de s a f ea si ble a n d eff e c2 tive app ro ac h to lif e p re dictio n a nd relia bilit y e sti matio n of p ro duct s w he n t he t e st ti me o r e x2 p e n s e a r e li m it e d.
K ey w ords : t y p e2Ⅰce n s o r e d dat a ; regre s sio n a n al y s i s ; e xt re me di s t ri b utio n ;
Wei b u ll di s t r i b utio n ; no r mal di s t ri b utio n ; locatio n2scale di s t r i b utio n ;
relia b ilit y e s ti matio n
回归分析在工程中有着广泛的应用,但传统的回归分析方法只适用于来自正态分布的完全数据,对于工程中常见的截尾数据情况却无法处理. 然而在高可靠性、长寿命产品的寿命试验中,由于受试验时间、经费等限制, 试验常常采取截尾处理,得到的均为截尾数据. 为此,文献[ 1 ]建立了定数截尾数据线性回归分析方法,将传统的回归分析推广到定数截尾数据的情况. 文献[ 2 ]提出了单
收稿日期:2009209230 ; 修订日期:2009210215
基金项目:国家自然科学基金( 10472006) ; 北京市教育委员会共建项目( X K100060418)
作者简介:傅惠民( 1956 - ) , 男, 浙江遂昌人,“长江学者”特聘教授, 博士, 从事小样本信息技术、软校准技术、数据融合方法及可靠性研究.
143 第1 期傅惠民等:定时截尾数据回归分析方法
状态定时截尾数据最佳线性无偏估计方法. 本文
进一步建立了定时截尾数据回归分析方法,解决
了上述难题.
与第q
i
+ 1 个顺序统计量的协方差, 均可通过公
式计算或查表[ 4 ] 得到.
取式(5) 中Q 关于a , b 和σ的偏导, 并令它们
等于零,求解可得参数a , b 和σ的估计值为极值分布定时截尾数据回归分析
极值分布是工程中常见的一种分布. 设y 为
服从极值分布的随机变量, 其分布函数为
1
^a= y- - ^bx- - σ^u-( 7)
L 22 L 1y - L 12 L 2y
^b =(8)
L 11 L22 - L22 1
L11 L2y - L12 L1y F(y) = 1 - e xp - e xp ( 1)σ^= ( 9)
L 11 L22 - L22 1
式中μ( x) 为状态x 下的位置参数,σ> 0 为尺度
参数. 令y°= [ y - μ( x) ]/σ,则y°为标准极值分布
随机变量. 设μ( x) 与x 之间成线性关系,即
式中
q
i
+ 1
m
- 1
y = ∑∑
3
νik l(10)
y ik
μ( x) = a + b x
则随机变量y 可表示为
y = a + b x+ε, ε～EV ( 0,σ)
式中a , b 和σ为不依赖于x 的待定参数.
11 1 参数估计
( 2)
u- (11)
x-
( 3)(12)
q +1
m i
3 ∑∑νikl(13)
n=
i = 1 k, l = 1
设在状态x i( i = 1 ,2 ,, m) 下有n
i 个试样进q i +1
m
∑∑νikl-
( x i-- ) ()( )行定时截尾试验, 到y 3时刻截止, 共有q ( 1 ≤L1y=x y i k-y14
i i
i = 1 k, l =1
q i < n i ) 个试件失效, 得到定时截尾数据y i1 ≤q +1
m i
y i2 ≤≤y i q, 它们可以看作来自极值分布大小∑∑νikl ( u ik--
u) ( y y )(15)
L 2y=--
il
i
i = 1 k , l = 1
为n i 的样本的前q
i
个顺序统计量Y i1 ≤Y i2≤
q
i
+1
m
≤Y i q 的一个取值. 根据文献[ 223 ] 可知,∑∑νikl ( x i- x- ) 2 (16)
L11=
i
= y 3可以看作该样本的第q 个区间统计i =1k, l = 1
y i ( q + 1)i i
i
q i +1
m
量Y i( q + 1) 的一个取值.
∑∑νik l ( x i-- )
x- ) ( u ik i( )
17
L12=-u 由式( 3) 可知
y i k = a + b x i+εik ,
i = 1 k , l =1
q i +1
m
εik ～EV (0 ,σ)∑∑νik l ( u ik- ) (- )
L22=-u u il-u( )
18
i = 1 k , l = 1
i = 1 , 2 ,
根据文献[ 223 ] , 令
m
, m ; k = 1 , 2,, q
i
+ 1(4)
可以证明, ^a,^b 和σ^是参数a , b 和σ的最佳
无偏整体估计. 理论及Mo n t e Ca rlo 模拟表明,当
m
q i +1
Q = ∑∑( y i k- a -bx i -
数据总量q = ∑q i + 1 稍大时,可认为它们近
i = 1 k, l = 1
i = 1
似服从正态分布.
11 2 位置参数与尺度参数的置信限估计
位置参数μ( x) = a + b x的点估计为
μ^(x) = ^a+ ^bx
其置信度为γ的单侧置信上、下限分别为
μup ( x) = ^a+ ^bx +
σu i k )νikl ( y i l bx i - σu i l )(5)
- a -
式中
- 1
[νikl ] (q + 1) ×(q + 1)= [νikl ] (q + 1) ×( q + 1)(6)
i i i i
且u ik ( k= 1 , 2 ,, q
i
) 为标准极值分布大小为n i
( 19)的样本的第k 个顺序统计量的均值,νik l ( k , l = 1 ,
, q
i
) 为标准极值分布大小为n i 的样本的第k
2 ,
个与第l 个顺序统计量的协方差, u i( q + 1) 为标准
i
极值分布大小为n i + 1 的样本的第q i + 1 个顺序σ^uγ
ω( 20)
uγ(c13 + c23 x) +
2
1 - uγc33
统计量的均值,νik ( q + 1)=νi ( q + 1) k( k = 1 , 2 ,, q
i
+
i i
1) 为标准极值分布大小为n i + 1 的样本的第k 个
y - μ( x)
σ
144
航 空 动 力 学 报 第 25 卷
式中
μlow ( x ) = ^a + ^bx +
σ
^u γ 2
1 ω ω P = + c 33 l n l n +
ω
( 21)
u γ ( c 13 + c 23 x ) - 1 - u 2 c 33
1 - P
γ 1
其中
2 c 1
3 + c 23 x l n l n
( 31)
1 - P
ω= u 2
( c 13 + c 23 x ) 2
+ γ 区间[ y P low , y P up ]为概率为 置信度为 2γ- 1 的置信区间 , 即
P 的百分位值 y P
( 1 - c 33 u 2 γ) ( c 11 + 2 c 12 x + c 22 x 2 ) ( )
22 尺度参数σ置信度为γ 的单侧置信上 、下限
分别为
P{ y P low ≤y P ≤y P up } = 2γ- 1
( 32)
σ^ Weibull 分布定时截尾数据回归
分析
设随 机变 量 t 服从 Wei b ull 分 布 , 分 布 函
数为
σup =
2 ( 23) 1 - u γ
σ^
c 33
σlow =
( 24)
1 + u γ
c 33
其中
m
( c j j ) 3 ×
3 = F ( ) = 1 - e xp t - , t > 0
1 2 η( x )
( 33)
式中 m > 0 为形状参数 ,η( x ) > 0 为特征寿命. 令
y = l n t , 则 y 服从极值分布 , 分布函数见式
( 1) , 且
μ( x ) = l n η( x ) q i +1
q i +1
q i +1
- 1
m
m
m
∑∑
νikl
∑∑
νikl
x i
∑∑
νikl
u ik
i = 1 k , l = 1 i = 1 k , l = 1 i = 1 k , l = 1 q i +1
q i +1
q i +1
m
m
m
∑∑
νikl
x i u
ik
i
i = 1 k , l = 1 i = 1 k , l = 1 i = 1 k , l = 1 q i +1
q i +1
q i +1
( 34)
m
m
m
1 σ=
m
∑∑
νikl
u ik
∑∑
νikl
x i u ik
∑∑
νikl
u ik
u il
i = 1 k , l = 1
i = 1 k , l = 1
i = 1 k , l = 1
此时 , Wei b u ll 分布定时截尾数据回归分析转化 为极值分布定时截尾数据回归分析问题进行 处理.
( 25)
区间 [μlo w ( x ) ,μup ( x ) ] 和 [σlo w ,σup ] 分别为
μ( x ) 和σ置信度为 2γ- 1 的置信区间 , 即
正态分布定时截尾数据回归分析
设随机变量 y 服从正态分布 , 分布函数为
3 P{μlow ( x) ≤μ( x) ≤μup ( x) } = 2
γ- 1
P{σlow ≤σ≤σup } = 2
γ- 1 ( 26)
( 27)
11 3 百分位值点估计及其置信限估计
随机 变 量 y 概 率 为 P 的 百 分 位 值 F ( y) =Φ ( 35)
y P = 式中μ( x ) 为状态 x 下的均值 ,σ为标准差 ,Φ( ·) 表示标准正态分布函数.
当μ( x ) 与 x 之间成线性关系时 , 有
1
μ( x ) +σl n l n
的无偏估计为 1 - P
1 ^y P = ^a + ^bx +σ^l n l n 1 - ( 28)
P y = a + b x +ε, ε～N ( 0 ,σ
2
) ( 36)
其置信度为γ的单侧置信上、下限分别为
根据 11 1 节可得到参数 a , b 和σ的最佳无偏整体
σ^
ν 估计 , 由式 ( 7) ～ ( 9) 给出 , 此时 , u ik ,
ikl 分别表示 标准正态分布下相应统计量的均值和协方差 , 可
y P up = ^a + ^bx + · 1 - u 2
c 33 γ 1
+ u 2
( c 13 + c 23 x ω l n l n
) + u γ [ 5 26 ]
通过公式计算或查表得到 .
均值μ( x ) 和标准差σ的点估计及置信度为γ 的置信限按 11 2 节方法计算. 随机变量 y 概率为 P 的百分位值的点估计及置信度为γ的置信限按
γ P
1 - P
( 29)
σ^ y P lo w = ^a + ^bx + ·
1 - u 2
c 33 γ 11 3 节方 法计算 , 由式 28 ～ 30 给 出 , 式 中
( ) ( ) 1 + u 2 ( c 13 + c 23 x ω ) l n l n - u γ γ 1 - P
P 1 l n l n 1 - 由Φ ( P ) 代替.
- 1 ( 30)
P
y - μ( x )
σ
t
145 第1 期傅惠民等:定时截尾数据回归分析方法
表1 M onte Carlo 模拟数据
T a ble 1 Da ta of M onte C arlo simulation 一般位置2尺度分布定时截
尾
数据回归分析
4
序号
i
温度
T i / K
样本量截尾时刻定时截尾数据
t i k / h 设随机变量
数为
y 服从位置2尺度分布, 分布
函
t 3 / h
n i i
1 343 6 1 000 883
2 36
3 6 500 221 , 490
- ∞< y < + ∞
F,(37)
3 383 6 300 72 , 232 , 246 , 286 式中μ( x) 为状态x 下的位置参数,σ> 0 为尺度
参数, F (·) 表示位置2尺度分布函数.
当μ( x) 与x 之间成线性关系时, 根据11 1 节
可得到待定参数的最佳无偏整体估计,由式( 7) ～
( 9) 给出,此时, u ik ,νikl 分别表示标准位置2尺度分
布下相应统计量的均值和协方差.
位置参数μ( x )和尺度参数σ的点估计及置
信度为γ的置信限按11 2 节方法计算. 随机变量
y 概率为P 的百分位值的点估计及置信度为γ的
置信限按11 3 节方法计算, 由式(28) ～(30) 给出,
4 403 6 100 29
令y = l n t ,则y 服从极值分布EV μ( T) ,σ,
且μ( T) = l nη( T) ,σ= 1/ m .采用本文第1 节方法对
其进行回归分析, 得到参数 a , b 和σ的最佳无
偏整体估计分别为
^a= - 81 828 6 , ^b = 5 6081 813 0 , σ^= 01 465 1 .
因此,形状参数μ( T) 的无偏估计为
μ( T) = l nη( T) = - 81 828 6 + 5 6081 813 0/ T
其置信度为2γ- 1 = 01 9 的区间估计列于表2 .
为了便于比较,本文还将表1 中的定时截尾
数据简单看成定数截尾数据, 采用定数截尾数据
回归分析方法进行分析, 所得的μ( T i ) 的置信区
间估计结果也列于表2 .
由表2 对比结果可以看出,两种方法所得的
μ( T i ) 的置信区间都包含其真值, 但本文方法所
得的各个温度下的置信区间长度都明显较短, 也
就是说,本文方法计算精度明显较高. 大量Mo n t e
Ca r lo 模拟证明, 在相同的样本量下, 本文方法所
得的置信区间长度与采用定数截尾数据回归分析
所得结果相比, 平均可以缩短20 %～40 % , 有效
提高了定时截尾数据的回归分析精度.
式中l n l
5 算
51 1
设某产品寿命t 服从Wei b ull 分布, 形状参
数m = 2 , 特征寿命η( T) 与绝对温度T 有如下
关系:
l nη( T) = a + b/ T
式中a = - 10 , b = 6 000 . 现采用Mo nt e Ca rlo 模拟
方法在4 个温度下对产品进行定时截尾试验,具
体模拟条件和模拟数据列于表1 .
μ( T i ) 置信区间估计的对比结果( 2γ- 1 = 01 9)
表2
Comparison results o f conf i d ence interv al s of μ( T i ) ( 2γ- 1 = 01 9)
T a ble 2
区间长度相对差值
本文方法定数截尾数据回归分析方法
μ( T i)真值 d2 - d1
×100 %置信区间区间长度d
1
置信区间区间长度d2d2
μ( T1 )
μ( T2 )
μ( T3 )
μ( T4 )
71 493
61 529
51 666
41 888
[ 71 027 9 , 81 386 6 ]
[ 61 334 6 , 71 205 5 ]
[ 51 549 0 , 61 312 5 ]
[ 41 657 4 ,51 692 1 ]
11 358 7
01 870 9
01 763 5
11 034 7
[ 71 099 8 , 91 202 0]
[ 61 278 4 , 71 526 8]
[ 51 290 7 , 61 278 7 ]
[ 41 061 3 ,51 494 2 ]
21 102 2
11 248 4
01 988 0
11 432 9
351 37
301 24
221 72
271 79
特征寿命为η( T) 的Wei b ull 分布,大量工程经验
表明η( T) 与绝对温度T 有如下关系:
l nη( T) = a + b/ T
51 2 工程算例
表3 给出某产品在4 个不同温度水平下的定
时截尾寿命试验数据[7 ] , 它们服从形状参数为m 、
y - μ( x)
σ
146
航 空 动 力 学 报 第 25 卷
表 3
定时截尾寿命试验数据
因此有
η
^ ( T ) = e xp ( - 22 . 215 7 + 11 593/ T )
m
^ = 1/σ^ = 1 . 170 5 温度 T 下位置参数μ( T ) 和尺度参数σ置信 度为γ的置信上、下限分别为 T a b le 3 T ype 2 Ⅰ c ensored l i f e t est d ata
序号 i
T 1 = 358 K 寿命/ h T 2 = 398 K 寿命/ h T 3 = 423 K 寿命/ h T 4 = 448 K
寿命/ h 1 465 150 41 2 2 776 167 44 3 1 854 3 u γ
0 μup ( x) = - 22
. 215 7 + 11 593/ T + ·
1 - 01 011 1 u 2
γ 3 1 546 179 48 7 ω
u γ (
- 0 . 045 8 + 19 . 240 4/ T ) + 4
5 6 7 2 490
3 000
4 241
5 415 192 250 251 270 51 58 69 80 7 8
10 12 0
1 854 3 u γ
μlow ( x ) = - 22
. 215 7 + 11 593/ T + ·
1 - 01 011 1 u 2
γ ω
u γ ( - 0 . 045 8 + 19 . 240 4/ T ) -
01 854 3 σup = 8 6 107 305 83 13 u γ 01 011 1
01 854 3
1 - 9 7 273 329 101 14 σlow =
10 > 19 930 341 105 16 1 + u γ 01 011 1
11 > 19 930 372 112 18 其中
12 > 19 930 391 129 21 ω= u 2 2 γ ( - 0 . 045 8 + 19 . 240 4/ T ) +
13 > 19 930 560 159 31 ( 1 - 01 011 1 u 2 2
γ) ( 2 . 699 1 - 2 221/ T + 459 343/ T )
14 > 19 930 630 161 34 温度 T 下置信度为γ、概率为 P 的对数寿命 y P 的点估计及置信上限 、下限分别为
15 > 19 930 702 176 35 16 > 19 930 707 184 40 1
^y P = - 22 . 215 7 + 11 593/ T + 01 854 3 l n l n
1 - P
17 > 19 930 754 185 42 01 854 3
18 > 19 930 823 190 43 ·
= - 22 . 215 7 + 11 593/ T +
y 1 - 01 011 1 u 2
P up γ 19 > 19 930 835 195 48 1 19 . 240 4 2
ω P
ln l n 1 - P + u γ ( - 0. 045 8 +
) + u γ 20 > 19 930 > 1 097 > 245 50 T
21 > 19 930 > 1 097 > 245 56 01 854 3 y P low = - 22. 215 7 + 11 593/ T + 1 - 01 011 1 u 2 ·
22 > 19 930 > 1 097 > 245 75 γ
ln l n 1 + u 2 19 . 240 4 23 > 19 930 > 1 097 > 245 > 90 ω P γ ( - 0. 045 8 + ) - u γ 1 - P T
24 25
> 19 930
> 19 930
> 1 097
> 1 097
> 245
> 245
> 90
> 90
式中
2
1 ω =ω+ 01 011 1 l n l n
+
P 1 - P 将表 3 中寿命试验数据进行对数变换 , 则变 换后数据服从位置参数为μ( T ) = l n η( T ) 、尺度 参数为σ= 1/ m 的极值分布. 采用本文第 1 节方 法对其进行回归分析 , 得到参数 a , b 和σ的最佳 无偏整体估计分别为
^a = - 22 . 215 7 , ^b = 11 593 , σ^
= 01 854 3 . 其协方差矩阵为
( - 0 . 091 6 + 38 . 480 8/ T ) l n l n 1
1 - P
结 论
6 1) 提出一种定时截尾数据回归分析方法 , 该
方法可以充分利用截尾时刻的试验信息 , 并将各 个状态下的定时截尾数据作为一个整体进行统计
分析 , 全面开发利用各个状态下的试验信息 , 因此 具有信息量大 、精度高的特点.
2) 详细讨论了工程中常见的极值分布定时
^a 11 970 0 - 8101 563 8
- 01 033 5
- 8101 563 8 - 01 033 5 ^b C o v = 335 261 141 043 1
141 043 1 01 008 1
σ
147 第1 期傅惠民等:定时截尾数据回归分析方法
截尾数据回归分析问题,给出回归系数和标准差的最佳无偏整体估计,以及百分位值的点估计及其置信限估计.
3) 建立适用于Wei b ull 分布、正态分布以及一般位置2尺度分布的定时截尾数据回归分析方法.
4) 本文方法对个别状态下仅有一个失效数据的情况也能很好地进行处理.
cal St rengt h ,2009 ,31 ( 6) :90529091 (i n Chi ne se)
傅惠民. 区间统计量及其分布[ J ] . 机械强度,2005 ,27 ( 6) :
75227571
FU Hui mi n . Int er val st ati stic a nd it s di st ri butio n [J ] . J o ur2
nal of Mecha nical St rengt h ,2005 ,27 ( 6) : 75227571 (i n Chi2
ne se)
Bala kri shna n N , Cha n P S. Or der st ati stics f ro m ext r e me
val ue di st ri butio n : Ⅰ- Ta ble of mea ns , va ria n ce s a nd co2
va ria nce s[J ] . Co mmunicatio n s i n St ati stics2S i m ul atio n and
Co mp ut atio n ,1992 ,21 ( 4) :1199212171
Pa r ri sh R S. Co mp uti ng exp ect ed val ues of no r mal o r der
st ati stics [ J ] . Co mmunicatio n s i n St ati stics2Si mul atio n and
Co mp ut atio n ,1992 ,21 ( 1) :572701
Par ri sh R S. Co mp uti ng va riance s a nd co va ria n ce s of no r2
mal o r der st ati stics[J ] . Co mmunicat io ns i n St ati s t ics2Si mu2
lat io n and Co mp ut atio n ,1992 ,21 ( 1) : 7121011
茆诗松, 王玲玲. 加速寿命试验[ M ] . 北京: 科学出版
社,1997 .
MAO Shi so ng , WA N G Li ngli ng .Accel erat ed lif e t e sti ng
[ M ] . Beiji ng : Science Pre ss ,1997 . (i n Chi ne se)
[ 3 ]
[ 4 ]
参考文献:
[ 1 ] 黄伟,傅惠民. 截尾数据线性回归分析方法[ J ] . 航空动力学报,2003 ,18 ( 4) :46524691
H U A N G Wei , FU H u i mi n . Li n ea r regressio n a nal y si s
met ho d f o r censo red dat a [J ] . J o ur nal of Aero space Po wer ,
2003 ,18 ( 4) :46524691 (i n Chi ne se)
[ 2 ] 傅惠民, 岳晓蕊. 定时截尾数据最佳线性无偏估计方法[J ] . 机械强度,2009 ,31 ( 6) : 90529091
FU Hui m i n , Y U E Xiao r ui . Be s t li nea r unbia sed e sti mat e
met ho d f o r t ype2 Ⅰcen so red dat a [ J ] . J o u r n al of Mechani2 [ 5 ] [ 6 ] [ 7 ]。