最近我在公司里做了一个小项目，帮助产品部门建立一个数据模型来预测产品的维修率和返修成本，其中有一步需要估计二参数威布尔分布的参数。在网上看了一些论文，威布尔参数估计的方法有很多种，比如常见的有极大似然估计法，最大相关系数优化法，最小二乘法等等。 因为考虑我目前仅是做模型验证，因此我采用了计算量相对比较小的最小二乘法进行估算，并在EXCEL里利用趋势图直接获得相关参数。
第一部分 威布尔分布（Weibull Distribution）介绍
威布尔分布是可靠性分析和寿命检验的理论基础。
威布尔分布是三参数的连续性的概率分布，其概率密度为：
f
(
t
,
m
,
η
,
γ
)
=
(
η
/
m
)
(
t
−
γ
η
)
(
m
−
1
)
e
−
(
t
−
γ
η
)
m
（
t
>
0
,
m
>
0
）
f(t, m, \eta, \gamma) = (\eta/m)(\frac{t-\gamma} {\eta})^{(m-1)}e^{-(\frac{t-\gamma}{\eta})^m} （t>0, m> 0）
f(t,m,η,γ)=(η/m)(ηt−γ)(m−1)e−(ηt−γ)m（t>0,m>0）
t 是自由变量，一般表示产品工作的时间，
η
\eta
η 是比例参数（scale parameter），
m
m
m是形状参数(shape parameter)，
γ
\gamma
γ是位移参数（location parameter）。
此处，我们研究产品的返修时间是以月为单位计时，从产品出厂第一个月就可能发生返修，因此我们可以假设
γ
\gamma
γ为0，概率分布退化为二参数形式：
f
(
t
,
m
,
η
)
=
(
η
/
m
)
(
t
η
)
(
m
−
1
)
e
−
(
t
η
)
m
（
m
>
=
0
）
f(t, m, \eta) = (\eta/m)(\frac{t}{\eta})^{(m-1)}e^{-(\frac{t}{\eta})^m} （m>= 0）
f(t,m,η)=(η/m)(ηt)(m−1)e−(ηt)m（m>=0）
当
m
>
1
m>1
m>1时， 图像如下图中红线，呈现为先递，增达到一个峰值后再递减。当
m
<
1
m<1
m<1时， 图像则呈现如下图中蓝线，呈现为单调递减函数。
我们在研究在保修期内的产品从出厂到返修中间间隔的月份数的分布时1，发现其曲线和威布尔分布当
m
>
1
m>1
m>1时的形态非常相似。
第二部分 最小二乘法拟合威布尔分布推导
我使用产品返修记录作为实验数据，来估算威布尔分布的两个参数 η \eta η和 m m m，首先要考虑的是威布尔分布是一个非线性的模型，我们想使用最小二乘数法的话，首先就要对它做一个线性化的变化，使之呈现 Y = A X + B Y=AX+B Y=AX+B的形式。
让我们先来看一下威布尔分布的累计分布函数：
F
(
t
)
=
1
−
e
−
(
i
/
η
)
m
(
t
>
0
,
m
>
0
)
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
(
1
)
F(t) = 1-e^{-(i/\eta)^m} (t >0, m > 0) ..............................(1)
F(t)=1−e−(i/η)m(t>0,m>0)..............................(1)
设有n个实验样本，按失效时间先后得到寿命数据
t
i
t_i
ti，对应的第
i
i
i个样本累计失效概率
F
(
t
i
)
F(t_i)
F(ti)，可以使用中位秩近似算法
F
(
t
i
)
=
i
−
0.3
n
+
0.4
F(t_i)=\frac {i-0.3}{n+0.4}
F(ti)=n+0.4i−0.3来估算。
在得到失效时间和累计失效概率估值后，即可使用最小二乘法来进行参数估计了。
首先，对(1) 式左右两边分别求两次对数得到：
l
n
(
l
n
(
1
1
−
F
(
t
)
)
)
=
m
l
n
(
t
)
−
m
l
n
(
η
)
ln(ln(\frac1{1-F(t)}) )= mln(t)- mln(\eta)
ln(ln(1−F(t)1))=mln(t)−mln(η)
令
y
i
=
l
n
(
l
n
(
1
1
−
F
(
t
i
)
)
)
y_i =ln(ln(\frac1{1-F(t_i)}) )
yi=ln(ln(1−F(ti)1))，
x
i
=
l
n
(
t
i
)
x_i=ln(t_i)
xi=ln(ti)，
a
=
m
a = m
a=m，
b
=
m
l
n
(
η
)
b=mln(\eta)
b=mln(η)， 原式可以变换成
y
i
=
a
x
i
+
b
y_i=ax_i+b
yi=axi+b。
我们需要找到这样一组（a，b）使得
L
=
Σ
i
=
1
n
(
y
i
′
−
y
i
)
2
L=\Sigma_{i=1}^n{(y_i'-y_i)^2}
L=Σi=1n(yi′−yi)2最小化，
y
i
′
y_i'
yi′是使用中位秩近似算法得到的
F
(
t
i
)
F(t_i)
F(ti)代入
l
n
(
l
n
(
1
1
−
F
(
t
i
)
)
)
ln(ln(\frac1{1-F(t_i)}) )
ln(ln(1−F(ti)1))得到的结果，的此处的
L
L
L可以看出来是方差和。
L
=
Σ
i
=
1
n
(
y
i
′
−
y
i
)
2
L=\Sigma_{i=1}^n{(y_i'-y_i)^2}
L=Σi=1n(yi′−yi)2
=
Σ
i
=
1
n
(
y
i
′
−
a
x
i
−
b
)
2
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
.
(
2
)
=\Sigma_{i=1}^n{(y_i'-ax_i-b)^2}........................(2)
=Σi=1n(yi′−axi−b)2........................(2)
(2) 是一个二次函数，通过求极值的方法可以求得其极值点。为了方便运算，我们采用矩阵表达方式:
A
=
(
a
b
)
A = \left ( \begin{matrix}a\\ b\end{matrix} \right )
A=(ab)
X
=
(
x
1
.
.
.
1
x
n
.
.
.
1
)
X = \left ( \begin{matrix} x_1...1\\ x_n...1\\ \end{matrix} \right )
X=(x1...1xn...1)
Y
=
(
y
1
.
.
.
y
n
)
Y=\left ( \begin{matrix} y_1\\ ... \\ y_n\\ \end{matrix} \right )
Y=
y1...yn
那么
(
2
)
(2)
(2)就可以表示成
L
=
(
Y
−
X
A
)
T
(
Y
−
X
A
)
L =(Y-XA)^T(Y-XA)
L=(Y−XA)T(Y−XA)，为求
L
L
L的最小值，对
A
A
A求导，得到
d
L
d
A
=
X
T
(
Y
−
X
A
)
=
0
\frac{dL}{dA}=X^T(Y-XA)=0
dAdL=XT(Y−XA)=0，求得
A
=
(
X
T
X
)
−
1
X
T
Y
A=(X^TX)^{-1}X^TY
A=(XTX)−1XTY。
由此，我们可以使用EXCEL来求解A。求得
A
=
(
a
b
)
A= \left ( \begin{matrix}a\\ b\end{matrix} \right )
A=(ab) 后，
m
=
a
,
η
=
e
b
m
m=a, \eta=e^{\frac{b}{m}}
m=a,η=emb 。
第三部分 使用EXCEL来计算
(1) A7:A56记录了样本序号，一共50个样本，分别标记为1…50。
(2) B7:B56将失效时间 t t t依次填入。
(3) C7:C56在C7填入公式 =LN(B7), 并扩展至C56。
(3) D7:D56全部填入1，和C7:C56构成了上一章提到的 X X X。
(4) E7:E56在E7填入公式 =C7*C7，并扩展至E56。
(5) F7:F56填写中位秩近似值，在F7中填入公式 =(A7-0.3)/(MAX($A$7:$A$56)+0.4) 并扩展至 F56。
(6) G7:G56填写Y值，在G7中填入公式 =LN(LN(1/(1-F7)))，并扩展至 G56。
(7) H7:I8是计算得到的 X T X X^TX XTX， H7的公式为 =SUMPRODUCT(C7:C56,C7:C56)， I7的公式为 =SUMPRODUCT(C7:C56,D7:D56)， H8的公式为 =SUMPRODUCT(C7:C56,D7:D56)， I8的公式为 =SUMPRODUCT(D7:D56,D7:D56) 。
(8) J7:K8是 ( X T X ) − 1 (X^TX)^{-1} (XTX)−1，使用数组函数 {=MINVERSE(H7:I8)} 得到2x2的矩阵。
(9) L7:M56是 ( X T X ) − 1 X T (X^TX)^{-1}X^T (XTX)−1XT，L7输入公式 =SUM(C7*$J$7,D7*$K$7) 并且扩展至L56。M7 输入公式 =SUM(C7*$J$8,D7*$K$8) 并且扩展至M56。
(10) N7是 a a a，输入公式 =SUMPRODUCT(G7:G56,L7:L56)。
(11) O7是 b b b，输入公式 =SUMPRODUCT(G7:G56,M7:M56)。
(12) B2是 m m m，输入公式 =N7。
(13) B3是 η \eta η，输入公式 =EXP(-O7/B2)。
(14) Q7:Q56用于做你和评估，Q列为输入
t
t
t，R列为
f
(
t
)
f(t)
f(t)，S列为样本数据的分布，作图显示
f
(
t
)
f(t)
f(t)和样本数据的匹配程度。
参考文献
[1] 史景钊 花恒明 李祥付 威布尔分布参数估计在EXCEL中的实现方法研究 中国科技论文在线 (http://www.paper.edu.cn) TB114.3
[2] 包小庆 刘志强 吴永忠 李冬梅 （2017）双参数威布尔分布函数的确定及曲线拟合 能源与环境 2017.NO.4 8-9
完
出厂当月记为0， 次月为1， 以此类推。统计每一台返修机的返修月份和出场月份的差获得间隔月份。 ↩︎
1266
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
