基于似然比检验的双截尾威布尔分布区间估计.pdf

第5 3 卷第6 期 2 0 17 年3 月 机械工程学报 J O U R N A LO FM E C H A N I C A LE N G I N E E R I N G V 0 1 5 3N O 6 M a r 2 0 17 D o I 1 0 3 9 0 1 册埴E 2 0 1 7 0 6 2 0 3 基于似然比检验的双截尾威布尔分布区间估计枣 南东雷贾志新李威 北京科技大学机械工程学院北京1 0 0 0 8 3 摘要 针对现有的双截尾威布尔分布模型估计过程中点估计无解析解和缺少参数区间估计方法的问题 提出基于极值求解思 想和蒙特卡罗方法的参数点估计数值求解方法和基于似然比检验理论的双截尾威布尔分布模型的参数区间估计方法 在计算 过程中 提出随机数取值区间多步细分的变区间迭代数值求解方法 结合某数控外圆磨床故障间隔时间数据使用双截尾威布 尔分布进行分析 并与两参数和三参数威布尔分布分析结果进行对比 实例结果表明 建立的变区间迭代算法和基于似然比 检验理论的区间估计方法正确可行 具有较好的估计效果 关键词 双截尾威布尔分布 似然比 参数估计 蒙特卡罗 中图分类号 T H l7 P a r a m e t r icI n t e r v a lE s t im a t io nf o rD o u b ly T r u n c a t e dW e ib u ll D is t r ib u t io nB a s e do nL ik e lih o o d R a t ioT e s t N A N D o n g le i J I AZ h ix inL IW e i S c h o o lo f M e c h a n ic a lE n g in e e r in g U n iv e r s it yo f S c ie n c ea n dT e c h n o lo g yB e ij in g B e ij in g1 0 0 0 8 3 A b s t r a c t S in c en oa n a ly t ic a ls o lu t io nw a sf o u n dint h ep r o c e s so f p a r a m e t r icp o in te s t im a t io na sw e lla ss o lu t io nf o rt h ep a r a m e t r ic in t e r v a le s t im a t io nf o rD o u b ly T r u n c a t e dW e ib u lld is t r ib u t io n a lle x t r e m ev a lu es o lu t io nm e t h o da n dM o n t eC a r lom e t h o db a s e d p a r a m e t r icp o in te s t im a t io ns o lu t io nm e t h o da n dlik e lih o o dr a t iot e s tt h e o r yb a s e dp a r a m e t r icin t e r v a le s t im a t io nm e t h o df o r D o u b ly T r u n c a t e dW e ib u lld is t r ib u t io nm o d e la r ep r o p o s e d D u r in gc a lc u la t io n t h em u lt i s t e pd e t a ile din t e r v a lc o n v e r s io nit e r a t iv e n u m e r ic a ls o lu t io nm e t h o df o rr a n d o md a t ain t e r v a lisp u tf o r w a r df u r t h e r O nt h eb a s iso fC N Ce x t e r n a lc y lin d r ic a lg r in d in gf a ilu r e d a t a D o u b ly T r u n c a t e d W e ib u lld is t r ib u t io nisu s e df o ra n a ly s isa n dc o m p a r is o nw it ht w o p a r a m e t e ra n dt h r e e p a r a m e t e rW e ib u ll d is t r ib u t io n A c c o r d in gt or e s u lt s t h ein t e r v a lc o n v e r s io nit e r a t iv ea lg o r it h ms e ta n dlik e lih o o dr a t iot e s tt h e o r yb a s e dp a r a m e t r ic in t e r v a le s t im a t io nm e t h o da r ec o r r e c ta n da v a ila b lea sar e s u lto f g o o de s t im a t io ne f f e c t K e yw o r d s d o u b ly t n m c a t e dW e ib u lld is t r ib u t io n lik e lih o o dr a t io p a r a m e t e r se s t im a t io n M o n t eC a r lo 0 前言 目前威布尔分布模型在可靠性分析过程中最 为常见 为计算方便多设位置参数为零将三参数威 布尔分布简化为两参数威布尔分布 继而使用最小 二乘法或极大似然估计参数点估计算法进行模型拟 合 并使用K S 假设检验或卡方检验对结果进行 检验 贾祥等 J 针对极少失效数据的参数估计问题利 用威布尔分布结合轴承失效数据对当前应用的分布 曲线法进行了改进 张根宝等u 1 使用改进的两参数 国家自然科学基金资助项I j 5 1 2 7 5 0 3 5 2 0 1 6 0 9 2 0 收到初稿 2 0 1 6 1 2 0 5 收到修改稿 威布尔分布模型对加工中心失效数据进行了分析处 理 并结合自助法完成对小样本数据的参数点估计 文献 3 利用两参数威布尔分布对材料疲劳可靠性 进行了分析并提出两种回归分析方法 文献 4 提出 了威布尔分布和对数正态分布的多线性回归算法 贾志新等 5 J 针对D K 7 7 线切割机床小样本问题提出 了基于K 均值聚类算法和径向基函数神经网络的扩 充算法 判定故障样本数据服从两参数威布尔分布 D A S l6 l从可靠性设计角度使用两参数威布尔分布对 蜂窝制造系统进行了分析 赵永翔等 7 1 基于两参数 威布尔分布描述了随机滚动接触疲劳应力幅和寿命 的关系 并对某型号轴承钢进行了可靠性分析 在 机床部件及系统可靠性分析中 威布尔分布被广泛 应用与研究 8 9 1 机械工程学报 第5 3 卷第6 期 在实际可靠性分析过程中故障间隔时间数据 并非从零开始也不为无穷大 而是在一定时间区间 内 例如为评价机床在可用期内的可靠性设定评估 时间区间从而剔除老练期故障数据 或进行一定时 间的老练试验从而降低使用期故障率等 lm J 而两 参数和三参数威布尔分布无法有效反映这些情况 在一定程度上影响了机床可靠性评价的准确性 W I N G O 1 0 l提出了双截尾威布尔分布模型参数点估 计方法 弥补了威布尔分布对双截尾数据分析的缺 陷 Z H A N G 等u 叫假设左截尾时间为零右截尾时间 未知 给出了右截尾情况下模型的参数估计数值解 法 然而此模型参数点估计无法取得解析解 且缺乏 有效的参数区间估计方法 从而影响分析的完整性 常用的区间估计方法有F is h e r 信息矩阵法 l引 B a y e s 方法 刘等 但F is h e r 信息矩阵法对样本量要 求较高且计算量较大 B a y e s 方法需要历史试验数 据以确定验前分布 而似然比检验方法并不存在这 些问题 l引 因此 本文提出基于似然比检验理论的 双截尾威布尔分布的区间估计方法 并基于蒙特卡 罗方法对模型参数进行点估计和区间估计数值 求解 1双截尾威布尔分布及其参数估计 1 1 双截尾威布尔分布模型 双截尾威布尔分布是威布尔分布的重要模型 之一 其累积分布函数及概率密度函数为 1 0 F x e x p 古 e X p f 古矽 e x p e x p T p 仆 一圆 1 弘矿 则双截尾威布尔分布模型的概率密度函数转 化为 a p x 肛1e x p 一o t x p 八为2 面i下葡匆 e x p I 一口u J e x p l一口 J 假设容量为n 的样本j 阢l抛 以 服从双截尾 威布尔分布 则模型似然函数为 三 口 x II 厂 口 p ig 尊柰端 模型对数似然函数为 m L a z ln 筇 ln 一口 矽一 一 础 1 e X p 一口 护一 6 分别对式 6 各参数求偏导数 并令其为零 求 解得 以 昙一喜m 薯 舀 丁型 巳 鹫匕出巫一 妒一 1 1 l毛 f f ln L z 2 ll l T p T p 1 1 p n T r pI nL 南一善n 矽一 一 1n Tp Lp expE丽 fl可 TP L 8 1 一e x p 西 矿一上尸 I 一 式 7 和式 8 没有解析解 但可以利用数值解法求 解 定义极值函数 G 2 看葛一善 一 一 n T p 五 L 再p e x 丽pI d 汗r f 1 可 T p L p 9 1 一e X p 卜西 F 卢一 l 并将式 7 代入式 9 中 即可得到以P 为白变量 的函数G 以I G I e 为目标函数使用蒙特卡 罗方法 1 即可求得符合要求的岔 其中e 为目标误 差值 代入式 7 即可得到对应的舀 从而进一步 得到台 1 3 似然比检验区间估计 似然比检验由L a w le s s 提出 在可靠性分布模 型区间估计中广泛应用 14 1 一 一 上丝上矿 一 一 一 I L 厂一 厂 2 0 1 7 年3 月 南东雷等 基于似然比检验的双截尾威布尔分布区间估计 2 0 5 似然比检验理论认为待估参数目和 的似然函 数L O 与极大似然函数L 6 分 比值的自然对数 的负二倍渐近于置信度为1 一O 自由度为1 的卡方 分布Z 1 即 一2 1 1 1 l坐L O 塑 f 1 L L Of 1 J l z 1 1 LJ 式 1 0 可改写为 硼 f 1 三 否 夕 e X p f 牮1 1 1 z I 在同一组样本中模型参数的极大似然点估计 值一般为唯一值 当置信度卜仃不变时 式 1 1 右 侧部分为常量 设 c 工 蛳 e x p f 一掣1 1 2 则式 1 1 可转化为 m t O f 1 I n C 1 3 通常情况下 在一定区间内式 1 1 中一个参数 确定时 另一参数有两个根随之确赳1 4 从而形成 两个参数的取值区间 然而该式没有解析解 可以 使用数值解法求解 令 o o ln C m L O 0 1 4 以式 1 4 为目标函数 取合适的区间和参数 即可使用蒙特卡罗算法求解 l5 将所有符合条件式 1 4 的解 曰 f 1 代入平均故 障间隔时间 M e a nt im eb e t w e e nf a ilu r e M T B F 的定 义式 M T B F I t f t d t 1 5 即可得到一组M T B F 值 由蒙特卡罗方法原理 可以认为符合条件的随机值足够多时 该组数据的 最小值与最大值即可视为M T B F 的置信区间 1 4 假设检验与拟合优度检验 在数控机床可靠性分析中常用的模型假设检 验方法有卡方检验和K S 检验 由于卡方检验只适 用于大样本而且容易犯I I 型错误 即错误的接受了 原假设 K S 检验更为精细且适用于小样本情况 因此本文采用K S 假设检验 做原假设H o 拟合的 分布模型F 0 等于理论分布F n t 将R 力与只0 的 最大偏差D 作为理论分布与经验分布的差异度作为 K S 检验统计量p J D s u p lE t F f I 1 6 式中R 忙模型理论分布 由中位秩公式确定 给定显著性水平O 通过查K S 检验极限分布 表 可得指定显著性水平的临界值元 比较统计量 D 是否小于允 若小于则接受原假设 否则拒绝 为测试模型拟合优度 定量表示应用不同方法 对同一样本分析的精确程度 引入相对均方根误差 N o r m a liz e dr o o tm e a ns q u a r ee r r o r N R S M E 作为评 价指标 川 N R M S E 1 7 式中 户 f f 模型计算累积失效概率值 由式 1 7 可得 参数估计效果越好 曲线拟合 精度越高 N R M S E 值越小 2 参数估计数值求解 2 1 参数点估计数值求解 本文使用长期跟踪的某数控外圆磨床5 1 个整 机故障间隔时间数据进行分析 其中模型左截尾时 间L 0 1h 右截尾时间T 30 0 0 h 如表1 所示 表1 某数控外圆磨床故障间隔时间数据表 序号i故障间隔时间他序号i故障间隔时间咖 l2 22 73 0 3 6 22 32 83 3 3 2 34 62 94 0 2 5 42 2 53 04 2 0 0 52 3 O3 l4 3 2 4 62 9 93 24 4 0 9 74 4 83 34 6 8 0 86 5 03 44 8 2 5 97 1 53 54 9 0 6 1 07 5 33 65 0 8 3 1 l7 7 O3 75 1 8 5 1 21 0 8 13 85 1 9 3 1 3 H 2 23 96 8 4 9 1 4 1 1 5 84 07 7 8 5 1 5 1 1 7 04 18 0 1 0 1 61 1 8 84 28 3 5 0 1 71 4 9 24 38 4 2 3 1 81 6 1 04 49 0 8 5 1 91 6 7 44 511 3 6 1 2 01 8 0 64 611 5 7 O 2 1 2 0 7 04 711 8 4 5 2 22 3 6 94 812 4 8 9 2 32 4 5 3 4 9l3 3 7 7 2 4 2 7 3 05 0l4 5 0 8 2 52 7 3 05 l29 8 9 0 2 62 9 2 6 2 0 6 机械工程学报 第5 3 卷第6 期 以p 为自变量 使用蒙特卡罗方法生成1 0 0 个 服从均匀分布的随机数据 代入式 9 中 绘制 G f 1 I G l图 如图1 所示 S o 豁 旧 坦 晕 i 岁 一 彤状参数 3 图1G D 和I G 励I 曲线图 通过程序初步判定局部极小值点 并视为疑似 解 如图1 所示 f G f 1 0 在 0 1 区间内有两个 疑似解 为获得足够精度 本文采用变区间迭代算 法 如图2 所示 图2 变区间迭代算法流程图 1 分别记录所有符合条件的拐点位置三个数 据净k 1 X i 2 X i 3 其中i 1 2 n 设初始计算 K N Y g A 声 2 设新的区间为阻 声羽 定义为 轴7 掣雨m J 掣 1 8 式中 迭代次数 m 区间控制参数 m O 1 决定疑似解 计算区间 取值越小速度越快 极值点 被忽略的风险越高 本文取m O 5 区间内再次使用蒙特卡罗方法生成服从均匀 分布的k 个随机数据代入式 9 中 得到新的 值 3 判断IG f 1 I e 是否成立 如果成立则继续 下一步 否则转到步骤 1 4 判断结果是否为唯一值 否则逐步缩小e 值 转到步骤 1 若为唯一值则结束迭代并输出最 终计算结果 其中步骤 4 为收敛条件 由于存在 基墨lG f 1 l 0 对于任意极小I T I 直e 若式 7 和式 8 有唯一解 则IG f 1 I 必然收敛于一点 将表1 中数据使用上述算法代入式 7 和式 8 迭代求解得西 o 0 0 6 0 矽 o 8 3 7 5 代入式 3 中得 毋 4 4 8 7 1 由式 1 可得该机床的故障模型为 o 9 9 99 一e x p 一0 0 0 60 x m 8 51 F z 2 蒜厉 j 双截尾威布尔分布模型拟合结果如图3 所示 o 籁 陋 妊 隶 器 嘴 0 5 0 0 0 0 0l5 0 020 0 025 0 030 0 0 故障间隔时间 h 图3 双截尾威布尔分布拟合曲线图 使用两参数威布尔分布拟合的分布模型为 O 8 3 6 1 F x 1 一e X p 志J 使用三参数威布尔分布拟合的分布模型为 F X 1 e x p I 一 等厂1 7 对三种模型做K S 检验 计算各模型M T B F 值 并进行N R M S E 拟合优度检验 结果如表2 所示 表2 假设检验和拟合优度检验 可靠性模型D DM T B F hN R M S E 两参数威布尔分布O 0 4 54O 3 2 404 6 8 9 10 0 3 57 三参数威布尔分布O 0 6 310 4 5 095 2 8 1 5O 0 4 72 双截尾威布尔分布O 0 5 10O 3 6 404 6 8 3 4O 0 3 54 2 0 1 7 年3 月 南东雷等 基于似然比检验的双截尾威布尔分布区间估计 给定显著性水平1 7 0 0 5 通过查K S 极限分 布表 可得旯 1 3 6 由表2 可得三种分布均符合 z D 旯 因此该数控外圆磨床故障数据服从两参 数威布尔分布 三参数威布尔分布和双截尾威布尔 分布 在拟合优度方面 双截尾威布尔分布的M T B F 值更接近实际观测值 N R M S E 值在一定程度上优 于两参数威布尔分布和三参数威布尔分布 因此双 截尾威布尔分布模型能更好地对该机床进行可靠性 分析 2 2 参数区间估计数值求解 本文使用蒙特卡罗方法求解双截尾威布尔分 布的参数置信区间 根据第2 1 节中参数点估计结 果分别确定适当的参数初始取值区间 即可进行数 值求解 取置信度1 一仃 0 9 5 1 依据上一节中参数点估计结果确定初始取 值区间为 O 6 5 1 p 3 0 0 6 5 0 2 确定参数 和0 的最佳区间 以确定声下 限为例 在0 的取值区间内使用蒙特卡罗方法随机 生成服从均匀分布的数据p 9 1 晓 岛 朗 和声 下限风一并代入式 1 4 若存在任一值可使 c o 0 成立 则扩大口下限 否则继续下一步 其中 和0 的上下限交替进行 直至确定卢和 0 的最终取值区间 计算过程如图4 所示 为清楚 显示效果 图4 中将G O 0 且lG O l e 的 0 隐藏处理 q 籁 诋 葵 漤 尺度参数0 图4 双截尾威布尔分布参数区间估计图 3 为节约计算时间和提高计算精度 取符合 条件G O 0 的区域上下左右四个极值点 分别 利用第2 1 节中变区间迭代算法求解 经过多次迭代计算 最终随机点取值区间分别 为 0 6 0 1 i0 0 3 0 0 7 0 0 在此区间内随机 生成1 0 5 个服从均匀分布数据 以式 1 4 为条件函 数 通过式 15 计算双截尾威布尔分布模型M T B F 的置信区间 计算结果如表3 及图5 所示 为清楚显示效果 图4 中将G O 0 且lG O I e 的 臼 值隐藏 处理 表3 极值点 尺度参数 图5 区间估计结果 由表3 可得 双截尾威布尔分布的尺度参数0 形状参数 和平均故障间隔时间M T B F 的区间估计 结果如表4 所示 表4 区间估计结果 模型参数参数的置信区间 尺度参数0 形状参数卢 平均故障间隔时间M T B F h 3 1 3 1 5 6 8 3 9 7 0 6 4 04 1 0 4 50 3 5 2 5 7 8 9 2 9 9 对两参数和三参数威布尔分布的形状参数卢和 平均故障间隔时间M T B F 进行区间估计分析 结果 如表5 所示 表5 不同模型区间估计结果对比 可靠性模型 参数置信区间 形状参数口平均故障间隔时间M T B F h 此数控外圆磨床故障数据的故障模型参数区 间估计结果显示 三参数威布尔分布的形状参数置 信区间精度最低 双截尾威布尔分布的形状参数置 信区间略低于两参数威布尔分布 M T B F 区间估计结果显示三参数威布尔分布 置信区间精度最低 双截尾威布尔分布精度低于两 参数威布尔分布 2 0 8 机械工程学报第5 3 卷第6 期 3结论 1 提出基于蒙特卡罗方法的双截尾威布尔分 布参数点估计数值求解算法 为提高运算效率及精 度引入变区间迭代算法 可以在确定初始区间后多 步细分以提高计算效率 2 提出基于似然比检验理论的双截尾威布尔 分布参数区间估计算法 并基于蒙特卡罗方法进行 数值求解 在确定整体取值区间后 根据极值点分 别使用变区间迭代算法求解 弥补了双截尾威布尔 分布无置信区间方法的不足 提高了使用该模型进 行可靠性分析的完整性 但估计精度有待进一步研 究改进提高 3 在此数控外圆磨床故障数据分析中 结果 显示双截尾威布尔分布的参数点估计精度较好 而 参数区间估计两参数威布尔分布的精度较高 使用 该算法和模型对某型号的中走丝线切割机床可靠性 分析也取得良好效果 参考文献 1 贾祥 王小林 郭波 极少失效数据和无失效数据的可 靠性评估 J 机械工程学报 2 0 1 6 5 2 2 1 8 2 1 8 8 J I AX ia n g W A N GX ia o lin G U OB o R e lia b ilit y a s s e s s m e n tf o rv e r yf e wf a ilu r ed a t aa n dz e r o f a ilu r e d a t a J J o u r n a lo f M e c h a n ic a lE n g in e e r in g 2 0 1 6 5 2 2 1 8 2 1 8 8 2 张根保 郭书恒 基于竞争威布尔模型的加工中心可靠 性评估 J 计算机集成制造系统 2 0 1 5 2 1 1 1 8 0 1 8 6 Z H A N GG e n b a o G U OS h u h e n g R e lia b ilit ye v a lu a t io n o fm a c h in in gt o o lc e n t e rb a s e do nc o m p e t in gw e ib u ll m o d e l J C o m p u t e rI n t e g r a t e dM a n u f a c t u r in gS y s t e m 2 0 1 5 2 1 1 1 8 0 1 8 6 3 B U T I K O F E RL S T A 恰R C Z Y KB R O O SM T w o r e g r e s s io n m e t h o d sf o re s t im a t io no fat w o p a r a m e t e r W e ib u lld is t r ib u t io nf o rr e lia b ilit yo fd e n t a l m a t e r ia ls J D e n t a lM a t e r ia ls 2 0 1 5 3 1 2 3 3 5 0 4 P I N A M O N A R R E ZM 0 R T I Z Y A N E ZJ w e ib u na n d lo g n o r m a lt a g u c h ia n a ly s isu s in gm u lt ip le lin e a r r e g r e s s io n J R e lia b ilit yE n g in e e r in g S y s t e mS a f e t y 2 0 1 5 1 4 4 2 4 4 2 5 3 5 贾志新 张宏斌 郗安民 基于径向基函数神经网络的 电火花线切割机床可靠性数据模拟生成 J 机械工程 学报 2 0 1 0 4 6 2 1 4 5 1 4 9 J I AZ h ix in Z H A N GH o n g b in X IA n n t in S im u la t in ga n d e x t e n d in gw ir ee le c t r ic a ld is c h a r g em a c h in in gr e lia b ilit y d a t ab yr a d ia lb a s isf u n c t io nn e u r a ln e t w o r k J J o u r n a lo f M e c h a n ic a lE n g in e e r in g 2 0 1 0 4 6 2 1 4 5 1 4 9 6 D A SK Ac o m p a r a t iv es t u d yo fe x p o n e n t ia ld is t r ib u t io n V Sw e ib u lld is t r ib u t io ninm a c h in er e lia b ilit ya n a ly s isma C M Sd e s ig n J C o m p u t e r s I n d u s t r ia lE n g in e e r in g 2 0 0 8 5 4 1 1 2 3 3 7 赵永翔 梁红琴 基于两参数W e ib u ll分布的概率疲劳 S N 曲线模型 J 机械工程学报 2 0 1 5 5 1 2 0 2 0 8 2 1 2 Z H A OY o n g x ia n g L I A N GH o n g q in M o d e lin go ft h e p r o b a b ilis t icf a t ig u eS Nc u r v e su s in gt h et w op a r a m e t e r w e ib u ll d is t r ib u t io n J J o u r n a lo fM e c h a n ic a l E n g in e e r in g 2 0 1 5 5 1 2 0 2 0 8 2 1 2 8 Y A N GZ K A NY C H E NF e ta 1 B a y e s ia nr e lia b ilit y m o d e lin ga n da s s e s s m e n ts o lu t io nf o rN Cm a c h in et o o ls u n d e r s m a ll s a m p led a t a J C h in e s eJ o u r n a lo f M e c h a n ic a lE n g in e e r in g 2 0 t 5 2 8 6 1 2 2 9 1 2 3 9 9 杨兆军 杨川贵 陈菲 等 基于P S O 算法和S V R 模 型的加工中心可靠性模型参数估计 J 吉林大学学报 2 0 1 5 4 5 3 8 2 9 8 3 6 Y A N GZ h a o j u n Y A N GC h u a n g u i C H E NF e i e ta 1 P a r a m e t e re s t im a t io no fr e lia b ilit ym o d e lo fm a c h in in g c e n t e rb a s e do np a r t ic les w a r lT Io p t im iz a t io na n ds u p p o r t v e c t o rr e g r e s s io n J J o u r n a lo fJ ilinU n iv e r s it y2 0 1 5 4 5 3 8 2 9 8 3 6 1 0 W I N G OD P a r a m e t r icp o in te s t im a t io nf o rad o u b ly 1 2 1 3 1 4 1 5 t r u n c a t e dw e ib u ll d is t r ib u t io n J M ic r o e le c t r o n ic s R e lia b ilit y 1 9 8 8 2 8 4 6 1 3 6 1 7 Z H A N GT X I EM O nt h eu p p e rt r u n c a t e dw e ib u ll d is t r ib u