Weibull分布拟合的经验似然方法

似然估计。
11 第三军医大学军事预防医学院卫生统计学教研室 ( 400038) 21 西南政法大学刑侦系 ( 400031 )
# 8#
C h inese Jou rnal ofH ealth Statistics , Feb 2009 , Vo. l 26, No . 1
E xp Exp
i 2 i
中国卫生统计 2009年 2月第 26 卷第 1期
# 7#
W eibull分布拟合的经验似然方法
郭波涛
1
王文昌
1
李辉智
2
刘
岭
1
易
东
1
=提 要 > 目的 利用经验似然方法拟合两参数 W e ibull分布。 方法 采用两 个矩条件 的经验似然 方法 , 同时利用 遗传算法进行求 解 , 应用于若干模拟数据集并与最大似然和分位数估 计的结果做比较。 结果 在大 样本时经验 似然估计 与最大似然估计 效果相当。 结论 经验似然方法可应用于 W eibull分布的拟合。 =关键词 >
- 2 /A - 1 /A
令 g 1i = x i - # ( 1 + 1 / A) K 2 / A) K G=
- 2 /A n n
, g 2i = x i - # ( 1 +
2 n j i= 1
, 构造相应的 Lagrange 函数:
E
log( npi ) - G
i= 1
Ep
i= 1
i
- 1-
E nr E p g
1. 3018 0. 8354 0. 5954 1. 0520 1. 8413
* : 最大似然法采用的是 M at lab71 0提供的 wb lfit函数计算 ; * * : 采用上下四分位数点进行估计
中国卫生统计 2009年 2月第 26 卷第 1期
# 9#
由表 1 可见基于大样本模拟数据 , 经验似然 + GA 计算的结果都是比较满意的, 而且同最大似然法估计 相似。分位数法为非参数估计方法 , 其估计值显然没 有最大似然和经验似然的效果好。 对于样本 1 而言, 最佳一次 GA 计算过程在第 50 个循环附近收敛 , 并且就种群进化而言也存在明显的 收敛趋势, 如图 1 所示。其他模拟样本也有类似的结 果。 进一步考察 , 通过经验似然估计的每个样本点的 概率值 p i 同实际 W e ib ull分布的累积分布的关系, 我 们同样利用样本 1 为例进行说明, 估计的累积分布同 理论分布之间的关系见图 2 , 可见两者有高度的相似 性。
表 1 大样本模 拟数据集三种方法的比较
真实值 样本号 1 2 3 4 5 A 1 0. 5 1 2 1 K 1 1 0 . 5 1 2 A 0 . 9906 0 . 5627 0 . 9968 2 . 1827 1 . 1281 最大似然 * K 1 . 1234 0 . 8788 0 . 5479 1 . 0745 1 . 7888 分位数法 * A 1. 3018 0. 5346 0. 9444 1. 9219 1. 1273
L( F ) =
F F ( {X } )
i i= 1
其中 F ( {X i }) 是分布 F 在 X i 处的概率质量 , 其中 i = 1 , 2 , , , n。 已知 X 1, X 2, , , X n 的经验累积分布函数 F n = n
- 1
E
n
i= 1 X i
I 使上式达到极大 , 即为 F 的非参数极大
n d 1 22
E
n i= 1
pi D 注意到 X 。 i
L(F n ) = n , 因而计算 R( H) 就成为关于 pi > 0 , i= 1 , 2 , , , n 求极大, 即: R( H) =
E
n p = 1, T (F p ) = H i= 1 i= 1 i
sup
F
n
( npi ) 。
同时 Owen
图 2 样本 1 计算累积分布与理 论分布比较
讨
论
本文以两参数的 W eibu ll分布 为例, 讨论 了如何 将经验似然应用于数据分布拟合, 同时进行了参数估 计。根据经验似然函数的特点, 提出了将遗传算法应 用于经验似然函数最优化的方案。通过对模拟数据的 验证, 可见经验似然估计的数值与一般的最大似然估 计相近。究其本源, 经验似然函数比实际上是一种截 面非参数似然比函数 , 在存在 W eibull分布 r 阶矩条件 约束条件下 , 其估计就赋予了参数化的信息, 由此以提 高其估计的准确性。可以推测, 样本量越大其估计的 值与似然函数越接近 , 估计的准确性越高, 反之, 若样 本量过少其估计效果不佳, 这一点我们也通过模拟数 据集做了验证, 具体结果见表 2 。
1 62 n
f =
E
i= 1
log 1 +
2
E
j= 1
+ rjg1 2
n
EP
m
m
Eg
i= 1
mi
/ 1+
E
j= 1
rj gji + P 3
E
1 n 1+
i= 1
2
- 1 rj gji
E
j= 1
其中 P i ( i = 1 , 2 , 3 ) 均为惩罚系数。 数据和结果 我们首先选取模拟样本量为 100 的样本数据, 分 别选择不同的 A和 K 值, 对每个样本分别进行 50 次 遗传运算, 其中初始种群大小为 50 , 最好种群数为 10 , 杂交率和变异率分别为 01 8 和 012, 选 择其中目标函 数最小的一次作为最终输出, 同时利用最大似然和分 位数估计, 结果见表 1 。
n 2 2
j= 1
s. t .
n
Eg
i= 1
1i
/ 1+
2
E
j= 1
, rj g ji = 0 = 0
E
i= 1
g 2i / 1 +
E rg
j= 1
j ji
对以上优化问题, M a tla b 采用序列二次规划法或 N ew ton 进行迭代求解 , 但这些方法对初始值的选择要 求比较高, 若初始值在真实值附近则能较快地收敛 , 反 之 , 计算时间过长或只能收敛到局部最优点, 这势必要 寻找对初值选择要求低 , 并保证搜索到全局或近似全 局最优点的算法。本研究中, 我们采用遗传算法解决 计算问题。 31 求解经验似然比函数的遗传算法构造 遗传算法
112
Ow en
定义如下经验似然比检验统计量 : L(F ) R( H) = sup | T (F ) = H ,FI F F L(F n )
1 32
其中, T 是分布 F 的某泛函 , 参数 H 是总体分布的特征。 由于 R ( H ) 的定义, 只有在观察值点有正概率质 量的离散分布才有可能使似然比达到极大, 于是为求 R ( H) 只需要限制在这样的离散分布类中求极大, 设分 布 F 在观察值 X i 有概率质量 pi > 0 , i= 1 , 2 , , , n, 并 记这样的分布函数 F 为 Fp , 即 F n =
1
j
E nr g
j= 1 2
9G 。 利用 E pi = 0得 G = n。 所以, 9pi i= 1
ji
n
1
。 最终 , 参数的求解转化为一个限
j ji
E rg
j= 1
制性最优化问题 : m in E log 1 +
iห้องสมุดไป่ตู้ 1 n 2
E
rj g ji
2
( 6) 关于初始种群 : 参数初始范围 的估计可以采 用分位点估计法得到 , 种群的大小要综合考虑发挥算 法的随机性 , 同时也要避免其盲目性 ; ( 7) 优化的目标函数 : 我们将约束 条件作为惩罚 项合并到目标函数中 , 得到如下函数 :
*
经验似然 + GA K A 0 . 9995 0 . 5857 1 . 0414 2 . 1435 1 . 1432 K 1 . 1305 0 . 9261 0 . 5818 1 . 0624 1 . 7878
GA 目标 函数值 0. 3115 1. 2127 1. 2586 0. 0422 0. 3127
1 32
证明在一定的条件下类似 W ilk s 定
理 R ( H) 服从卡方分布, 这样很容易将经验似然应用 于统计推断中。 21 W eibull分布参数估计的经验似然函数构造 两参数的 W e ib ull分布的密度函数形式为: f (x ) A - 1 A = KA x exp( - Kx ), 其中 K , A> 0 , 定义域为 x \ 0 。 W e ibu ll分布的形状取决于参数 A值, 因此, 称 A为 / 形 A 状 0 参数。 其生存函数形式为 S ( x) = exp ( - Kx ), 危 险函数为 h ( x) = KA x 其中 # ( A) =
i i
= # ( 1 + 1 / A) K
- 1 /A
,
2
经过若干代反复演化后 , 它将收敛到全局最优或近似 全局最优。遗传算法的优点是算法具有很好的稳健性 和全局最优性。 本研究利用 M atlab7 . 0 GA 工具包, 结合实际情 况采用如下的遗传算法来实现。
i ji 1 42
= # ( 1 + 2 / A) K
0 A - 1
法
, r阶矩为 [ # ( 1 + r /A) ] K ,
- r /A
u Q
]
A - 1 -u
e du 为伽玛函数。
利用经验似然方法进行拟合的关键在于如何加入 限制条件。 一个直观的想法在于 / 已知 0 每个样本的观 察值和其概率, 矩条件比较容易构造。 但对于参数分布 族如二项式分布、 正态分布等矩条件的参数不止一个, 如一阶为均值、 二阶为方差、 三阶为偏度等等。 若矩条 件的个数较多, 则参数估计限制过强 , 而少条件则不足 以估计参数。 在应用矩条件进行约束的同时也得到了 相应参数的经验似然估计, 可见经验似然的矩条件在 拟合中有两个作用: 限定数据的分布形式和估计相应 参数的值。 在本研究主要考 虑两参数的 W e ib ul l分布 情形, 由此选择 2 个矩条件 , 即:
W e ibull分布是 W. W eibull在 1939 年研究材料破 坏强度时建立的一种连续型分布, 因其有灵活的变换 形式, 简单的生存函数、 危险 函数和概率密度 函数形 式 , 因此广 泛应用于 可靠性分 析和生存 分析的 建模 1 12 中 。该分布是作为 连续分布中最小样本的 极限分 布而出现的 , 因此 , 在某些特定情况下 , 有人建议使用 W eibull分布作为极值的近似分布
-n
。
然而对于实际数据而言, 我们尚不清楚其分布形 式 , 或对其参数形式不能完全确定, 因此有必要采用一 些方法来拟合其分布 , 在本文的研究中我们尝试使用 经验 似然 ( Emp irica l L ik elihood , EL ) 来拟 合 W e ib ull 分布, 并将其结果同最大似然和分位数估计作一比较。 利用经验似然估计的求解除极少数分布外 ( 如分位数 估计的经验似然函数 ), 其他多参数分 布的经验似然 估计是很难直接求解的 , 甚至用一般的数值优化方法 1 22 都无能为力 。由于其自身的复杂性 , 当利用一般的 数值方法 (如 N etw on 法或拟 Ne t w on 法 ) 求它们的经 验似 然估计时 , 常常 导致不收 敛或收敛 到局部 极大 ( 极小 ) 值甚至鞍点 。所以本文拟采用遗传算法求 解经验似然估计值 , 将 GA 应用于经验似然的求解尚 未见报道。 方 11 经验似然 设 X 1, X 2, , , X n I é 独立有共同的累积分布 F, 则 F 的非参数似然是
是一种模拟生 物在自然环境 中的遗
传和进化过程而形成的一种自适应全局优化概率搜索 算法, 它是利用某种编码技术作用于成为染色体的编 码串 (可以是二进制串, 也可以是其他形式的码串, 如 实数编码、 区间值编码等 ), 并在选择、 杂交、 变异等算 子的作用下 , 按 / 适者生存 , 不适者淘汰 0的进化规则,
j= 1
其中 G , r1, r2 均为 Lagrange 系数。 优化问题转化为关于 {p i } 和 G , r1, r2 的 Lagrange 9 G 函数最大化。 固定 G , r1, r2, 对 pi 求导 , 并令 = 0 ,则 9pi 有 pi = G+ pi = n 1+
n 2
( 1) 编码 : 采用自然形式的实数编码, 每个染色体 包含两个基因, 分别表示参数 A和 K的一个估计; ( 2) 算术杂交算 子: M atlab 中 Crossover 的 cross2 overheuristic参数; ( 3 ) 非 均 匀 变 异 算 子: M atlab 中 M utat ion 的 G aussia n 参数 ; ( 4) 最好种群选择策略, 即从当前 父代和子代种 为选择适应值最大的若干个染色体作为新一 代染色 体; ( 5) 杂交率和变异率: 在作每一代的进化运算时, 按照杂交率和变异率来确定哪些染色体要做 杂交运 算 , 哪些基因需要作变异运算;