Gumbel分布的简单贝叶斯估计

2011年 3 月第 30 卷 第 2期 数理统计与管理 Jou rnal of A pp lied S tat isties an d M anagem en t M a r . 2 0 1 1 V 6 1. 30 N o .2 文章编号: 1002一1566 (2011)02一0234一06 G u m b el 分布的简单贝叶斯估计 鲁万波 焦鹏 (西南财经大学统计学院 , 四川 成都 61 0074) 摘要:在实际应用中, 两参数 G um bel 分布的贝叶斯估计往往需要预先知道 G um bel 参数的二维 联合先验分布. 由于获取先验分布的主观性和统计推断的复杂性 , 目前有关 G um bel 分布贝叶斯 估计理论及其性质的讨论还比较少 , 更不要说获得较为简单的 G um bel 分布的贝叶斯枯计. 本文 基于 K ajn ins kiy 和 V 舫iliy [l] 提出的简单贝叶斯估计过程 , 利用可靠度函数估计的区间形式示 先脸信息, 从而得到两个参数 G um be l分布的简单贝叶斯估计�基于此先脸信息, 该估计过程构 造了G um bel 参数的连续联合先验分布, 给出了在给定任意时.汽的可靠度 (或累积密度)及其标 准差的后脸估计 , 为可靠性与风险评估中简单快速的使用贝叶斯枯计,.1 画桩端辜件提供 了可能. 关锐词: 可靠性与风险评估; 贝叶斯估计;B et a 分布;G um bel 分布 中图分类号 : 0 212 .8 , 0 213 .2 文献标识码: A A S im P le B ay e sia n E st im at io n o f th e G u m b e l D istr ib u t io n L U M /a n一 b o J IA O P en g (Sch oolof Statisties, Southw estern U n ive rsity of F inanee and E eonom ies, S iehuan C hengd u 610074, C hina) A b straet: In Praetiee , it 15 neeessary to know a two-- d im ension aljoint Prior distribution of the G um - b el P a ra m eters in a d va n ee fo r m ak in g B ay esian estim at io n o f t~ p ajr am eter G u m b el d istrib u tion . It 15 relat iv ely fe w re su lts to d iseu ss th e th eo ry a n d n at u re o f B ay esia n estim 就 io n o f t一 P ar a m eter G u m b el d istrib u tio n , n o t to m en tio n g ain in g th e sim p ly B ay esian estim at ion . T h is P a P er o b ta ins B ay esia n e sti- m at io n of two-- P ar am ete r G u m b el d istrib u tio n b y u sin g a sim p le B ay esian estim at io n P ro eed u re p ro p o sed b y K am in sk iy a n d y 砧 iliy [1]. B as ed o n th is p rior in fo rm at ion , w h ieh ean b e p rese nted in th e fo rm o f the interval as sessm ent of the reliab ility fu netion , the Proeedure al low s eonstrueting the eontinuou s joint p rio r d istrib u tion o f G u m b el P aram eters as w ell as th e P o sterio r e stim at es of th e m ean an d stan d ar d deviation of the estim at ed reliability fu netion (or the eum ulative density fu netion) at any give n va lue of the exPosure va riab le. It 15 Possible to quiekly depiet extrem e eve nt by using the sim P le B ay esian estim 乱 ion in reliab ility a n d risk as sessm en t. K e y w o rd s : reliab ility an d risk as sessm en t , B ay esia n estim at io n , B eta d istrib u tio n , G u m b el d istrib u tion 0 引言 G um bel 分布 , 也称为 I型极值分布 , 是可靠性技术中一种常用的较复杂的分布类型 , 因 其对各种极端事件试验数据具有较好的拟合和描述能力 , G umb el 分布被广泛应用于可靠性 收稿 B 期: 2009 年 i 月 20 日 收到修改稿日期: 2009 年 12 月 22 日 荃金项 目:本文得到四川省统计科学研究计划项目 (2008sc 093) �西南财经大学 �211工程�三期青年教师成长 项 目 (21 IQ N 09020 )和西南财经大学 �,21 1工程三期 �统计学重点学科建设项 目资助. 鲁万波 , 焦鹏: G u m b el 分布的简单贝叶斯估计 2 3 5 和风险评估当中 {2一41�如今 , 贝叶斯估计技术已经被广泛的用于抽样检验 �可靠度分析以及 后验分布的估计 �但是 , 当我们考虑两个参数 G ullll�el分布的贝叶斯估计时 , 往往需要预先 知道 G um be l参数的二维联合先验分布 , 但是在实际运用中常常会遇到困难 , 如先验信息很 难得到或者统计推断过程过于复杂 �不利于应用 �一个 自然的想法是考虑使用共扼先验分布 , 这对于一些熟悉的常用已知分布 , 如正态分布或指数分布 , 是一个不错的思路 , 因为我们知 道 , 正态分布和指数分布都有各 自的共扼先验分布 , 且它们在实际应用中非常方便 �易于接 受. 与正态分布或指数分布不同 , 对于两个参数的 G um be l分布 , 当分布的两个参数都是随机 变量时 , 怎样给 G um be l分布选择合适的连续联合先验分布是十分困难的 � C ol eS an d P ow ell ls] 指出当利用贝叶斯估计方法进行统计推断的时候通常会遇到困惑 , 额外的先验信息有可能是 潜在的 , 仅仅给予某些极端行为已有的先验信息 , 这似乎不太合理 �理论上 , 一些学者讨论 了怎样建立 G um be l分布的贝叶斯估计 �Lye , H ap uar ac hch i和 R yan [0] 用两个熟知的贝叶斯 逼近过程: Li nd ley 和 Ti em ey & K ad an e 得到了 G um be l分布可靠度函数的贝叶斯估计 . Al i M ou sa , Ja he en 和 A hm adl 7} 给出了基于记录值的 G um be l分布的两个参数的贝叶斯估计过 程 �st 即hen so n 和 Jhw n[s] 把贝叶斯推断方法应用到了一个扩展模型上 , 这个模型可以把一个 非零概率精确地分配到 G um be l类型上 � 在可靠性和风险评估中怎样获得两个参数 G um be l分布的简单贝叶斯估计? K am insk iy 和 V白6ihy 提出 川, 先验信息可以通过可靠度函数估计的区间形式表示 , 并基于这个先验信息给 出了 W 匕ibul l分布的一个简单的贝叶斯估计过程. 在 K am insk iy 和 从访iliy 的启发下 , 我们考 虑了另外一种极值分布 一G um be l分布的简单贝叶斯估计过程 �这个过程利用可靠度函数估计 的区间形式表示先验信息, 构造 G unlbel参数的连续联合先验分布 , 并给出了在给定任意时刻 的可靠度 (或累积概率密度) 及其差的后验估计 �最后 , 我们举例说明了在可靠性与风险 评估中 , 怎样简明快速的运用该贝叶斯估计过程刻画极端事件 , 获取相应的贝叶斯估计结果 � 估计过程 假设被检验产品的寿命服从 G um be l分布 , 其概率密度函数和累积概率密度函数如下: f (t{a ,口)= F (艺}�,口) = 台二p(一宁 exp(一eXp(- 一 eX P t 一 a 口 艺一 a 口 )) (1) 一 oo < 云< + oo , (竹刀 其中 �(一oo < �< + oo )和 口, > 0) 分别为位置参数和尺度参数 � 不失一般性 , 假设已知两个预先给定的时点为 t:和 t: (亡: < 亡2), 以及在这两个时点的先 验信息一 累积概率密度函数的估计值歹蔽)和它的标准差, (可动), k一1,2�如图1所示,累积概率密度函数估计的不确定性由B et a 分布所描述. 该先验信息可以通过两种方法得到: 一种是基于大量数据的 K ap lan 一M ei er 和 G ree nw oo d 的非参数估计方法方法 , 该方法已经广泛使用;另一种方法是利用经验贝叶斯推断 , 此时先验信息通常由专家经验得到 日�基于这些先 验数据 , 我们可以得到一系列 G um be l分布的先验信息和后验信息 �在先验信息方面 , 我们可 以通过抽样的方式获得在任意固定时点 艺�的 G um be l参数的联合先验分布 �先验点估计和被 估计的累积概率密度函数的先验分布 �在后验信息方面 , 类似于先验信息的情况 , 我们同样可 以通过抽样的方式获得在时点 如的 G um be l参数的联合后验分布 �后验点估计和被估计的累 积概率密度函数的后验分布 , 其中 G ulllbel参数的后验点估计为联合分布最高后验概率密度 (众数), 类似于贝叶斯推断中的 H P D 区域. 2 3 6 数理统计与管理 第 30 卷 第 2 期 2011 年 3 月 刀乏F仓�),x0 (t:), , �(人)}} 厂司万吞万}/ 一尸洲万又动 豁曾圣戮亨 /减F 俩),凡(t:),, �伪)} r矛 勺 田 1 在给定时点 t�的先验信息 L l 先验分布 假设在固定时点获得的样本数据 , 要么是成功 , 要么是失败 , 为二元形式的观测数据 , 则 B et a 先验分布可以被用来描述被估计的累积概率密度函数的变化: �(;)-一华吐一 ,一�(1一,)几�一 �, 1 Lx o )I L几0 一 忿0 ) n �> 0, n �一xo > 0 , (2) 其中 , p 为随机变量 , 表示被估计的累积概率密度函数 F (t�), k 二1, 2;二0, n �一xo 是 B et a 分 布的形状参数 , 且 x �, n �为正;r( �)是通常的伽玛函数.基于任意给定时点 t�的先验累积概 率密度函数的估计值歹又万, 以及如图1中所示的标准差, (歹又)), 每个先验Bet �分布的参数可以通过矩估计法获得 110] : 一 尸 碑 户 . 气 ~ ~ 一 一 一 ~ ~ . ~ x �(艺�) = F (艺�) 一F (亡�), n �(艺�) x�(t*)= 二二二之二二F (t�) (3) 此外 , 每个先验B et a 分布的参数也可以通过极大似然估计法获得.在任意给定时点t�, n �(忿�) 基于先验累积概率密度函数估计值 F (t*), k = 1, 2 的对数似然函数为: 一�� ,一�凡�一��刃(淤�{��,万1�弩 �无),)一h刃(��({{{一塑({�),). 一 � (�)一气苏0气乙无) 一上j上11厂吓k) 十又no气不无)一x o又不无)一 1) In 气1 一岁杯k)) , 则未知参数二�(t*)和��(t �)的极大似然估计王苏动和而不又)是下面非线性方程组的解: 干王)��!{儿)!一生)��)��)�一���{�龄 (�一丽,,� 望L几o又不无)) 一 岁气x oL不儿)) = 一In 厂吓无), (5) 其中, (习一三镖华一箫 为digam ma函数�对刊溅性对数似然方程组(s) , 不存在解析解 , 可用数值方法求解.解非线性方程组的数值方法很多 , 如牛顿法 , 最速下降法和多元二分 法 [l1]等 �牛顿法具有收敛快 �形式简单的优点 , 但对初值要求较高 , 最速下降法计算简单 � 收敛性好 , 但收敛速度较慢 �戴家佳等 111]将解一元方程的二分法扩展到求解多元非线性方程 组 , 给出了多元二分法的算法 , 对 B et a 分布的参数极大似然估计进行模拟 , 发现当两个形状 参数的比增大时 , 无解的比例增加 , 随样本容量的减小而降低 �由于本文中可选用矩估计做为 极大似然估计的初始值 , 且在时点 t�的样本仅为单点样本 , 两参数估计植相差较大 , 我们使 用牛顿法获取非线性对数似然方程组 (5) 的解 . 鲁万波 , 焦鹏:G u lllb el 分布的简单贝叶斯估计 2 3 7 通过从被估计的 B et a分布中进行随机抽样 , 我们可以得到在给定时点 tl和 tZ 时 , 被估计 累积概率密度函数的n 对随机实现 {只(tl),只(亡2)},乞= 1, 2, � , n �对于任意一对{只(tl),凡(艺2)}, t: < tZ, 由于 F 川 是单调递增的 , 一定有 F (tl) < F 仕2), 它们各自的 G um be 分布的尺度参数 和位置参数可以通过下面公式得到: 亡2 一 艺i F *(tl)一F *(tZ) 艺IF *(tZ)一亡ZF *(tl)� 一万硕面丁不瓦不 (6) 这里, F *(��) = In (一In F (t�)), 无= 1,2, 并且 a = 一t�+ 口F *(��)�很明显, (6)式中n 对 G um be l分布参数 {a �,仄} 可以从 n 对随机模拟的 {凡(tl),凡(tZ)} 中得到 , 从而可以构造它们 的联合先验分布 �进一步 , G um be l累积概率密度函数 只(粼�,,风)的实现可以在给定的时点 如 外插获得 �利用 (3)式 , 这个分布可以再次用参数为 x�(t3)和 no(t3)的 B et a 先验分布逼近 � L Z 后验分布 在给定时点 亡�的 G um be l累积概率密度函数的后验分布: 在时点 t3 的 B et a 先验分布可 以通过贝叶斯估计获得 �与先验分布的情况类似 , 考虑 N 次试验有 : 次失败 , 每一次实验要 么是成功 , 要么是失败 , 因此在时点 艺3 可以得到二元形式的观测数据 �这样 , 后验分布也是 B et a 分布 , 其参数为: xo(亡3)* = x�(t3)+ :, no(t3)�= n �(艺3)+ N . (7) G umb el 分布参数的后验分布:给定时点 亡:, 亡:和 t3 (t3 是外插时的新时点), 从被估计的 B et a 分布中进行随机抽取 , 可以通过模拟的方式得到 G um be l参数的联合后验分布. 已知任 意固定的时点 t�, 无= 1, 2 和在 t3 时点累积概率密度函数的分布 , 由 (6)式可以得到先验累积 概率密度函数歹蔽)和它的标准误差, (歹蕊)), 每个Gumb el 累积概率密度函数的后验尺度 参数和位置参数也可以由 (6)式得到 �类似概率纸法 , 我们可以得到 G um be l分布参数的估计 值. G um be l参数的后验点估计就是联合分布最高后验概率密度 (众数)� 2 数值研究 假设一个机械零件的寿命服从 G um be l分布 , 专家估计在 tl二1, t: = 3 年零件发生故障 的可能性分别为成 百一0.01, 瓜 了一005 , 估计误差均为10% , 则发生故障概率标准差的估计分别是, (面又下)一0.001和, (歹获百))一�.005 .首先, 考虑Beta先验分布参数的估计�由(3) , 对给定时点 艺1 和 tZ, B et a 分布参数的矩估计分别为 {x � �(艺,{x � �(亡2 二 9 4 .9 5 n � �(亡2) 二 18 99 }; 以矩估计为初始值 , 由 = 98 �99 , n � �(tl) = 9899}) 式获得非线性方程组的 )临 , B eta 分布参数的极大似然估计分别为 {x: �(, 1) = 99.49,二: �(tl) = 9599.46}和 {x; �(tZ) = .45 ,7石不)一189 9.49 }�由于机械零件发生故障均为整数, 在给定的时点亡, 和,2, Bet a分布 解式和95 参数的估计分别是 {二�(艺1) = 99 , n �(tl) = 9899}和 {x�(艺2) = 95 ,二�(tZ) = 1899} �比较矩估计和 极大似然估计 , 前者计算简单 , 具有解析解 , 同时使用了零件发生故障的可能性和误差信息 ; 后者没有解析解 , 只有数值解 , 仅需要零件发生故障的可能性信息就可以估计 B et a 分布的参 数 ;在本文所考虑的问中 , 两者的估计较为接近 � 然后 , 从两个被估计的 B et a 分布中随机抽取累积概率密度函数 , 获得 n = 10000 个随机 实现 , 充分利用公式 (6), 可以得到 G umb el 参数的联合先验分布 , 如图2 所示 �与参数联合分 布的众数 (H PD )相对应的先验估计为{娜r= 5.879 ,岛二一4. 521} �若时点t3为第5年, 发生 2 3 8 数理统计与管理 第 30 卷 第 2 期 2011 年 3 月 故障的概率和它的标准差分别是 0.281 5 和 0.0388 �由公式 (3)或公式 (5) 可以得到在时点 如 的 Beta 分布的参数估计 {二�(亡3) = 37,n �(亡3)= 133}. 进一步 , 考虑在 t3 得到的数据 , 假设此时这些被观测的 12 个机械零件中 , 有 1 个出现 故障 �由公式 (7), 在 亡3 时点 , 后验 Bet a 分布参数可以观测到 , 它的估计值为 {亡�(云3)�= 38 加而亩沐一145 }�充分利用公式(6), 已知 �2和,3, 从被估计的Beta分布中抽取随机实现 � = 1000 0 , 就可以得到 G umb el 参数的联合后验分布 , 结果在图 3 中表示出来. G umb el 参数 的后验HpD 估计值为:{d;�s�一0.097,岛�:�= 2.513}� /// ����� HHHHH P OOOOO 图 2 G um bel 参数的联合先验分布 (10% 误差) 嘴乃 图 3 2 办 定名 盆0 盆6 确力 4 j 肠力 G um bel 参数的联合后验分布 (10 % 误差) 由于此时后验分布已经知道 , 给定 t3 , 故障概率 (累积概率密度函数) 的后验估计值和它 的标准差分别是 0.2620 和 0.0359 �数值检验和估计结果列在表 1 中 �本文的所有计算与作图 均使用统计软件 R 完成 � 表 1 先验与后验估计 (10% 误差) 先 先先验信息 息 后验信息 息 ttttti = l (给定))) 亡2 = 3 (给定 ))) 亡3 = 5 (推断 ))) t3 = 5 (推断 断 FFF (忿 ))) 0.01000 0.05000 0.281555 0 .262000 aaa (F (t)))) 0.00111 0.00555 0.038888 0.035999 ddddd 5 .8 7999 0 .09 777 声声 声 4.52111 2.81333 最后 , 为评价估计方法的优劣 , 我们对 G umb el 分布参数的联合先验分布和联合后验分 布进行了灵敏度分析 �还是用上面例子中的数据 , 如果把第 1年和第 3 年零件发生故障可能 性的误差增加到 20 % , 与预期的一样 , G um be l参数的联合先验和后验分布的二维散点图会比 10 % 误差的时候变得更宽; 同样 , 如果把第 1 年和第 3 年零件发生故障可能性的误差降低到 5% , G um bel 参数的联合先验和后验分布的二维散点图会比 10 % 误差的时候变得更窄 �限于 篇幅 , 我们在这里省略了相关的图示 , 如果读者需要 , 请与我们联系 �灵敏度分析的数值检验 和估计结果如表 2 所示 �对比表 1 和表 2 , 我们发现 , 无论是把误差增加到 20 % , 还是降低 到 5% , 在 t3 时刻的 F (t) 的先验信息和后验信息的推断结果变化并不大;而当把误差增加到 20 % 时, t3时刻的a( F (t) )的先验信息和后验信息的推断结果均增加了约两倍, 当把误差降低 到5% 时, 云3时刻的�(F (f) )的先验信息和后验信息的推断结果均减少了约两倍, 这与图示的 变化是一致的 �在 G ulllbel分布参数的先验估计与后验估计方面 , 位置参数 d 变化较大 , 尺度 鲁万波 , 焦鹏: G u m bel 分布的简单贝叶斯估计 2 3 9 参数应变化较小�以上灵敏度分析是在误差增加和减少两倍的情况下进行的, 如果误差的变 化不是特别明显 , 整体推断的结果会是比较稳健的 , 具有一定抵御外界干扰的能力 , 因此利用 这个方法进行估计 , 不仅简单 , 而且比较稳健 � 表 2 灵敏度分析 222220% 误差 差 5% 误差 差 先 先先验信息 息 后验信息息 先验信息 息 后验信息息 艺 艺艺i 二 111 亡 2 = 333 艺 3 = 555 t3 = 555 艺 i = 111 tZ = 333 t3 = 555 t3 = 555 (((((给定))) (给定))) �推断 ))) (推断 ))) (给定 ))) (给定 ))) (推断 ))) (推断 ))) FFF (t))) 0.01000 0.05000 0.283999 0.234 111 0.0 1000 0.05000 0.281444 0.276666 ��(F (亡 )))) 0.00222 0.01000 0.075444 0.060444 0.000555 0.002555 0.019555 0.018999 ddddd 6 .11222 0 .04 666 6 .12000 0 .0 1 111 口口 口 4.65999 2.70000 4.64000 2 76666 3 结论 本文考虑了 G um be l分布的一个简单贝叶斯估计过程 , 并进行了灵敏度分析 , 它的新颖之 处在于先验信息可以用可靠性函数估计的区间形式表现 , 而无须预先知道 G um be l参数的二 维联合先验分布 ;优良性是在可靠性与风险评估中, 相比于其它方法可以简单快速 �稳健的使 用贝叶斯估计刻画极端事件.该过程基于专家的先验信息 , 构造了G umb el 参数的连续联合先 验分布 , 并给出了在给定任意时刻的可靠度 (或累积概率密度)及其标准差的后验估计 �在理 论上 , 进一步的研究可以扩展到建立其它极值分布的简单贝叶斯估计 , 并且把它们统一到同 一个贝叶斯框架中 �在实际应用中 , 进一步的研究可以应用这个贝叶斯估计过程快速高效的 研究产品寿命的抽样检验计划 , 也可以应用于研究气象科学 �水文科学 �地震趋势预报等.后 续的相关研究正在进行中 � l参考文献 ] 11� M ark p K a而 nskiy and V曲 iliy V K rivtsov . A sim ple proeedure fo r Bayesian estim ation of the W 七 ibull distrib ution [J}.IE E E 肠 an sac tions on R eliab ility, 2005, 54(4): 612一616 � 冈 罗纯 ,王筑娟 �G um be l分布参数估计及在水位资料分析中应用 闭.应用概率统计, 2005 , 21 (2): 169 一17 5. [a] 陈云, 张志华�高可靠性软件分组数据的极值统计分析 15] .计算机应用与软件, 2006 , 23( 7):3于41 . [a] 冯平 ,王仲压, 田为民.基于二维 Gum bel 分布的长距离输水系统水文风险评估 [s] .灾害学, 2008 , 23 (l)= 23 es~26 . !5 � C oles 5 G an d p ow ell E A . B ayesian m eth od s in extrem e val ue m odelling: a review and new deve l- opm ents [J � .InternationalStat istieal Re view , 1996, 64(1): 119一136 �[6 � Lye M , H 即 uar ac hch i K p , and R yan 5. B ayes estim at ion of the extrem * va lue reliab ility fu netion [J{.IE E E 升 an saetions on Re liability, 1993, 42(4): 64 1芍44. !71 A liM ousa M A M , Jah een Z F , and A hm ad A A .B ayesian estim at ion, predietion and eharaeterization fo r the gu lnbel m od elb as ed on reeords �Jl.Stat isties, 2002, 36(1): 65一74. �8� Stephenson A and lh w n J .B 叮esian infe renee fo r extrem es: aeeount ing fo r the three extrem altypes �J �.E xtrem es, 2004, 7(4): 291一307. �9 � A yyub B M .E lieitat ion of E xpert O p inions fo r U neertai nty an d R isks � M l. B oea R aton, F L : C R C P ress , 20 0 2. �10} M odar ess M , K am ink siy M , an d K rivtsov V .Re liability E ngineering and R isk A nalysis: A p rac tieal G uide [M � .N ew Yo rk : M ar eelD ekker Ine., 1999.l川 戴家佳, 杨爱军, 杨振海. 极大似然估计算法研究 [s] .高校应用数学学报 (A 辑), 2009 , 24(3): 275一280 .