威布尔分布相关系数法参数估计在EXCEL中的实现方法研究

（2）
根据失效时间和累计失效概率即可用各种方法对其参数进行估计。
3. 最大相关系数优化法
对（1）式做变形处理，并取两次自然对数得到：
ln ln 1 = m ln(t − γ ) − lnη m 1− F (t)
（3）
令：
y = ln ln 1 , x = ln(t − γ ), B = lnη m
1− F (t)
n
∑( xi yi − nx ⋅ y) i =1
n
n
∑ ∑ ( xi2 − nx 2 )( yi2 − ny 2 )
i =1
i =1
（5）
式中：
∑ ∑ x
=
1 n
n i =1
xi
，
y
=
1 n
n i =1
yi
由(5)式可知，相关系数 R (x, y)是位置参数 γ 的函数，我们要求的 γ 应是使 R (x, y)最大
n
− ny 2 )
i =1
x − xi ti − γ
所以：
∑ ∑ ∑ ∑ ∑ ∑ n
uv '− vu ' = 2(
i =1
yi2
n
− ny 2 )(
i =1
xi yi
⎡n − nx ⋅ y) ⎢(
⎣ i=1
n
xi2 − nx 2 )
i =1
y − yi ti − γ
n
−(
i =1
xi yi
n
− nx ⋅ y)
i =1
i =1
i =1
则有：
dR(x, y)2 = uv'−vu'
dγ
u2
即需满足：
uv'−vu'= 0
因为：
∑ ∑ ∑ v ' =
dv dγ
=
n
2(
i =1
xi yi
− nx
⋅
y) d dγ
n
[
i =1
ln(ti
− γ ) yi
−
y
n i =1
ln(ti
− γ )]
∑ ∑ n
= 2(
i =1
xi yi
量)为 t1 ≤ t2 ≤ Λ ≤ tn ，对应的累计失效概率(经验分布函数)为 F (t1) ≤ F (t2 ) ≤ Λ ≤ F (tn ) 。 其中到第 i 个产品失效时的累计失效概率 F (ti ) 可用中位秩算法求得：
-1-
http://www.paper.edu.cn
F (ti )
≈
i − 0.3 n + 0.4
和$I$8>=0；
-4-
http://www.paper.edu.cn
⑶ 单击“求解”按钮，即可获得最大相关系数下的位置参数 γ =20.2395，如图 2 所示，
此时可获得最大相关系数 R(x,y)=0.99950878；
图 2 规划求解结果
3）使用图表功能求形状参数 m 和尺度参数 η ⑴ 插图散点图，横坐标为 xi，纵坐标为 yi； ⑵ 在散点图上添加趋势线，回归模型选择“线性”，并选择“显示公式”和“显示 R2 值”； ⑶ EXCEL 自动绘制回归直线，并把结果显示在图上，结果如图 3 所示。其中斜率 1.8486
ห้องสมุดไป่ตู้
n
−(
i =1
xi yi
n
− nx ⋅ y)
i =1
x − xi ti − γ
=0
（6）
满足（6）式的 γ 即为所求的位置参数。
4. 用 EXCEL 进行参数估计
（6）式所表示的方程十分复杂，解该方程一般是通过编程，用数值解法求出 γ ，然后
求再用最小二乘法或其他方法求解形状参数和尺度参数。MS EXCEL 具有强大的统计和计 算功能，其“规划求解”功能更是求解最优化问题的强有力工具，（6）式所表示的方程利用 EXCEL 的“规划求解”功能可很容易解出，然后再利用其散点图的趋势线功能即可求出形状 参数和尺度参数。本文通过实例，就相关系数优化法，用 EXCEL 进行求解。
于是（3）式化为线性方程
y=mx-B
（4）
令：
xi
=
ln(ti
−γ
); yi
=
ln ln 1−
1 F (ti )
,i
= 1,2,Λ
,n
则在 γ 已知时可利用最小二乘法进行回归求出 m 和 B， 进而求出尺度参数
η = exp( B ) 。 m
x 与 y 间的相关系数 R (x, y)为：
R(x, y) =
[2] Hallinan A J. A review of the Weibull Distribution[J]. Journal of Quality Technology. 1993,25(2) ：85-93 [3] Adatia A．，Chart L. K． Robust. Estimators of the 3-parameter Weibull Distribution[J]. IEEE
⑸ 在 D2 单元格中输入公式“=C2*C2”，用填充柄填充 D3～D6 单元格，D2～D6 单元
格的值为为 xi2 ；
⑹ 在 E2 单元格中输入公式“=(A2-0.3)/5.4”，用填充柄填充 E3～E6 单元格，E2～E6 单
元格的值为为 F (ti ) ，这里 F (ti ) 采用中位值算法，即 F (ti ) = (i − 0.3) (n + 0.4) ；
t 是产品的工作时间， t ≥ γ 。
当 m<1 时， 由式( 3 ) 给出的失效率是递减型的，适合于建模早期失效；当 m=1 时， 失效率为常数，即退化为指数分布，适合于建模随机失效；当 m>1 时，失效率是递增的， 适合于建模磨耗或老化失效。
设有 n 个产品进行寿命试验数据，按失效时间先后得到的寿命数据失效时间(顺序统计
n
∑ ⑽ 在 C8 单元格中输入公式“=SUM(D2:D6)”，C8 单元格的值为 xi2 ； i=`
⑾ 在 C9 单元格中输入公式“=AVERAGE(F2:F6)”，C9 单元格的值为 y ；
n
∑ ⑿ 在 C10 单元格中输入公式“=SUM(G2:G6)”，C10 单元格的值为 xi ⋅ yi ； i=`
例：选取 5 台某产品进行可靠性试验，失效时间分别是 27,32,36,42,49，已知产品寿命 服从威布尔分布，试估计分布参数。
1) 准备数据表 按图 1 准备数据表
图 1 位置参数估计数据准备表
⑴ 在 A2～A6 单元格中输入产品失效的顺序号 1～5； ⑵ 在 B2～B6 单元格中输入产品的失效时间 27、32、36、42、49；
∑ ⒂ 在 F8 单元格中输入公式“=SUM(H2:H6)”，F8 单元格的值为 n x − xi ；
i=1 ti − γ
∑n
⒃ 在 F10 单元格中输入公式“=SUM(I2:I6)”，F10 单元格的值为
y − yi ；
i=1 ti − γ
⒄ 在 I10 单元格中输入公式“=(C8-5*C7^2)*F10-(C10-5*C7*C9)*F8”，I10 单元格的值
（3）失效概率有中位秩算法、平均秩算法，采用的算法不一样，估计结果也会稍有不 同。
-5-
http://www.paper.edu.cn
参考文献
[1] Weibull W. A statistical distribution function of wide applicability[J] . Journal of Applied Microelectron Reliability. 1951,28(4)：613-617
i =1
x ti
− −
xi γ
⎤ ⎥ ⎦
n
n
∑ ∑ 因为 0<|R|≤1，所以 ( yi2 − ny 2 ) ≠ 0 ，否则 R=∞； ( xi yi − nx ⋅ y) ≠ 0 ，否则 R=0。
i =1
i =1
因此只有:
∑ ∑ ∑ ∑ n
(
i =1
xi2
n
− nx 2 )
i =1
y − yi ti − γ
2. 三参数威布尔分布模型
威布尔分布的寿命分布函数由下式给出
F (t )
=1−
⎡ exp⎢−
⎢⎣
⎜⎜⎝⎛
t
−γ η
⎟⎟⎠⎞m
⎤ ⎥ ⎥⎦
t ≥γ
(1)
式中：m 称为形状参数，m>0；η 称为尺度参数，η>0； γ 称为位置参数，也称最小寿
命，表示产品在 γ 以前不会失效，对于产品寿命有 γ ≥ 0 , γ =0 时退化为二参数威布尔分布；
1. 引言
威布尔分布是瑞典物理学家 Weibull W.分析材料强度时在实际经验的基础上推导出来 的分布形式[1]，国内外大量研究表明，用三参数威布尔分布比用对数正态分布往往能更准确 地描述结构疲劳寿命或腐蚀损伤的概率分布[2]，物理意义更加合理；在以损耗为特征的机械 零件寿命评估中，采用三参数威布尔分布比采用二参数威布尔分布拟合精度更高。因此，三 参数威布尔分布在强度与环境研究领域及机械零件磨损寿命评价中得到越来越广泛的应用。 在农业机械的强度设计中也经常要用到威布尔分布。
-3-
http://www.paper.edu.cn
⑶ 在 I8 单元格中输入位置参数 γ 的迭代初值，初值可选择接近于第一个失效时间，也
可用图估计法的估计值作为初值；
⑷ 在 C2 单元格中输入公式“=LN(B2-$I$8)”，用填充柄填充 C3～C6 单元格，C2～C6
单元格的值为 xi ，即 ln(ti − γ ) ；
就是（6）式的左边。
其他文字仅用于说明，与求解关系不大，可以不填。
2）使用“规划求解”功能估计位置参数 γ
⑴ 选择“工具→规划求解”功能打开规划求解参数对话框，目标单元格设为$I$10，目 标值设为 0，可变单元格设为$I$8:$B$2。
⑵ 对于产品寿命有 0 ≤ γ < t1 ，故单击约束条件“添加”按钮，添加约束条件：$I$8<=$B$2
即为形状参数 m，而 5.5088 即为 lnη m ，故η = exp⎜⎛ 5.5088 ⎟⎞ = 19.69 ⎝ 1.8486 ⎠
图 3 回归结果
5. 结语
（1）三参数威布尔分布的参数估计较二参数威布尔分布的参数估计更为复杂，用 MS EXCEL 解决三参数威布尔分布的参数估计问题，实用方便。
（2）用“规划求解”进行位置参数的估计时，要注意选择好迭代精度和初值，若精度和 初值选择不合适，可能得不到最满意的解，建议初值尽可能离第一个失效时间近一些，精度 不低于 10-7。
n
− nx ⋅ y)
i =1
y − yi ti − γ
-2-
http://www.paper.edu.cn
∑ ∑ ∑ u '
=
du dγ
=
n
(
i =1
yi2
−
ny 2 )
d dγ
n
[
i =1
ln2 (ti
−γ)−
1 n
⎛ ⎜⎝
n i =1
ln(ti
−
γ
)
⎞2 ⎟⎠
]
∑ ∑ n
= 2(
i =1
yi2
⒀ 在 H2 单元格中输入公式“=($C$7-C2)/(B2-$I$8)”，用填充柄填充 H3～H6 单元格，
H2～H6 单元格的值为 x − xi ； ti − γ
⒁ 在 I2 单元格中输入公式“=($C$9-F2)/(B2-$I$8)”，用填充柄填充 I3～I6 单元格，I2～
I6 单元格的值为 y − yi ； ti − γ
威布尔分布参数估计方法有很多， 国内外一直有人在进行相关研究[3-8]，现有几十种参 数估计方法，但多数只能用于形状参数和尺度参数的估计。在众多的估计方法中，能用于三 参数估计的并不多，见诸文献的有极大似然估计法、最大相关系数优化法、概率权重矩法、 灰色估计法、图估计法等，除图估计法外，其他方法大都计算复杂，应用不便，即便是计算 机水平发达的今天，也只能通过 Matlab 或其他计算机语言编程计算。EXCEL 提供了超强的 数学运算、统计分析等实用程序 ，利用它的规划求解功能可以快速、高效地求解三参数威 布尔分布的参数估计问题。
http://www.paper.edu.cn
威布尔分布参数估计在 EXCEL 中的实现方法研究
史景钊，花恒明，李祥付
河南农业大学机电工程学院，郑州（450002 ）
E-mail: jingzhaos@gmail.com
摘 要：三参数威布尔分布的参数估计比较复杂，大多数估计方法都需要编程计算。推导了 相关系数优化法求解三参数威布尔分布位置参数的公式，并把该公式利用 MS EXCEL 强大 的规划求解进行求解，实现了位置参数估计的非编程求解，同时利用图表功能实现了形状参 数和尺度参数的非编程求解。实例证明该法方便、快捷，可以满足估计精度。 关键词：可靠性；威布尔分布；参数估什；EXCEL 中图分类号：TB114.3
的那一个。根据求极值的方法，使 d R(x, y) / dγ = 0 的 γ 即为所求， d R(x, y) / dγ = 0 等
价于 dR2 (x, y) / dγ = 0 。
若令：
n
n
n
∑ ∑ ∑ v = ( xi yi − nx ⋅ y)2 ， u = ( xi2 − nx 2 )( yi2 − ny 2 )
⑺ 在 F2 单元格中输入公式“=LN(LN(1/(1-E2)))”，用填充柄填充 F3～F6 单元格，F2～
F6
单元格的值为为
yi
，即
ln
ln
1
−
1 F
(ti
)
；
⑻ 在 G2 单元格中输入公式“=C2*F2”，用填充柄填充 G3～G6 单元格，G2～G6 单元格
的值为 xi ⋅ yi ； ⑼ 在 C7 单元格中输入公式“=AVERAGE(C2:C6)”，C7 单元格的值为 x ；