韦伯分布(Weibull)参数矩估计MATLAB实现
二参数韦伯分布概率密度函数
f ( x ) = β η ( x η ) β − 1 e − ( x η ) β , β > 0 , η > 0 , x ≥ γ ≥ 0 f(x)=\frac{\beta}{\eta}\left(\frac{x}{\eta}\right)^{\beta-1} e^{-\left(\frac{x}{\eta}\right)^{\beta}}, \beta>0, \eta>0, x \geq \gamma \geq 0 f(x)=ηβ(ηx)β−1e−(ηx)β,β>0,η>0,x≥γ≥0
其中, β \beta β是形状参数, η \eta η是尺度参数。
β 和 η \beta和\eta β和η的矩估计
1. Weibull分布的k阶矩:
μ k = ( 1 η β ) − k β Γ ( 1 + k β ) \mu_{k}=\left(\frac{1}{\eta^{\beta}}\right)^{-\frac{k}{\beta}} \Gamma\left(1+\frac{k}{\beta}\right) μk=(ηβ1)−βkΓ(1+βk)
其中,
Γ
\Gamma
Γ是gamma函数:
Γ
(
s
)
=
∫
0
∞
x
s
−
1
e
−
x
d
x
,
(
s
>
0
)
\Gamma(s)=\int_{0}^{\infty} x^{s-1} e^{-x} d x,(s>0)
Γ(s)=∫0∞xs−1e−xdx,(s>0)
推导过程:
E
X
k
=
∫
0
∞
x
k
⋅
β
η
β
x
β
−
1
e
−
(
x
n
)
β
d
x
=
∫
0
∞
β
η
k
(
x
η
)
k
+
β
−
1
e
−
(
x
n
)
β
d
(
x
η
)
=
∫
0
∞
η
k
(
x
η
)
k
e
−
(
x
n
)
β
d
(
x
η
)
β
=
∫
0
∞
η
k
x
k
β
e
−
x
d
x
=
η
k
Γ
(
k
β
+
1
)
\begin{array}{l} E X^{k}=\int_{0}^{\infty} x^{k} \cdot \frac{\beta}{\eta^{\beta}} x^{\beta-1} e^{-\left(\frac{x}{n}\right)^{\beta}} d x \\ =\int_{0}^{\infty} \beta \eta^{k}\left(\frac{x}{\eta}\right)^{k+\beta-1} e^{-\left(\frac{x}{n}\right)^{\beta}} d\left(\frac{x}{\eta}\right) \\ =\int_{0}^{\infty} \eta^{k}\left(\frac{x}{\eta}\right)^{k} e^{-\left(\frac{x}{n}\right)^{\beta}} d\left(\frac{x}{\eta}\right)^{\beta} \\ =\int_{0}^{\infty} \eta^{k} x^{\frac{k}{\beta}} e^{-x} d x \\ =\eta^{k} \Gamma\left(\frac{k}{\beta}+1\right) \end{array}
EXk=∫0∞xk⋅ηββxβ−1e−(nx)βdx=∫0∞βηk(ηx)k+β−1e−(nx)βd(ηx)=∫0∞ηk(ηx)ke−(nx)βd(ηx)β=∫0∞ηkxβke−xdx=ηkΓ(βk+1)
2. 因此可以计算出一阶矩和二阶矩
m 1 = μ ^ k = ( 1 η ) 1 β Γ ( 1 + 1 β ) m 2 = μ ^ k 2 + σ ^ k 2 = [ 1 η ] 2 n { Γ ( 1 + 2 β ) − [ Γ ( 1 + 1 β ) ] 2 } \begin{aligned} m_{1} &=\hat{\mu}_{k}=\left(\frac{1}{\eta}\right)^{\frac{1}{\beta}} \Gamma\left(1+\frac{1}{\beta}\right) \\ m_{2}=\hat{\mu}_{k}^{2}+\hat{\sigma}_{k}^{2} &=\left[\frac{1}{\eta}\right]^{\frac{2}{n}}\left\{\Gamma\left(1+\frac{2}{\beta}\right)-\left[\Gamma\left(1+\frac{1}{\beta}\right)\right]^{2}\right\} \end{aligned} m1m2=μ^k2+σ^k2=μ^k=(η1)β1Γ(1+β1)=[η1]n2{Γ(1+β2)−[Γ(1+β1)]2}
让m2除以m1的平方,我们可以得到:
σ
^
k
2
μ
^
k
2
=
Γ
(
1
+
2
β
)
−
Γ
2
(
1
+
1
β
)
Γ
2
(
1
+
1
β
)
\frac{\hat{\sigma}_{k}^{2}}{\hat{\mu}_{k}^{2}}=\frac{\Gamma\left(1+\frac{2}{\beta}\right)-\Gamma^{2}\left(1+\frac{1}{\beta}\right)}{\Gamma^{2}\left(1+\frac{1}{\beta}\right)}
μ^k2σ^k2=Γ2(1+β1)Γ(1+β2)−Γ2(1+β1)
这正好是变异系数的平方,上式开方就是变异系数:
C
V
=
Γ
(
1
+
2
β
)
−
Γ
2
(
1
+
1
β
)
Γ
(
1
+
1
β
)
C V=\frac{\sqrt{\Gamma\left(1+\frac{2}{\beta}\right)-\Gamma^{2}\left(1+\frac{1}{\beta}\right)}}{\Gamma\left(1+\frac{1}{\beta}\right)}
CV=Γ(1+β1)Γ(1+β2)−Γ2(1+β1)
3. 我们可以借助变异系数CV来估计 β \beta β
- 首先我们定义一组候选 β \beta β,如 { β ∣ 0.1 < β < 10 } \{\beta|0.1<\beta<10\} {β∣0.1<β<10}的数组,根据这一系列候选 β \beta β计算出对应的 C V CV CV数组
- 然后根据数据的均值和方差计算数据的变异系数 C V ^ \hat{CV} CV^
- 找出 C V CV CV数组中和 C V ^ \hat{CV} CV^最接近的值,其对应的 β \beta β就是我们需要的估计值
4. 尺度参数 η \eta η可以利用下面的公式估计
η ^ = { x ˉ / Γ [ ( 1 / β ^ ) + 1 ] } \hat{\eta}=\{\bar{x} / \Gamma[(1 / \hat{\beta})+1]\} η^={xˉ/Γ[(1/β^)+1]}
其中
x
ˉ
\bar{x}
xˉ是数据的均值。
下方参考文献给出尺度参数的计算公式如下:
η
^
=
{
x
ˉ
/
Γ
[
(
1
/
β
^
)
+
1
]
}
β
˙
\hat{\eta}=\{\bar{x} / \Gamma[(1 / \hat{\beta})+1]\}^{\dot{\beta}}
η^={xˉ/Γ[(1/β^)+1]}β˙
这应该是一种笔误,在实际测试时按照上述公式计算的不对,翻阅其他文献(如这篇)是没有指数项的。
MATLAB实现
function [beta,eta]=estimateweibullMOM(vec)
beta_vec=0.01:0.001:10;
g1=gamma(1+2./beta_vec);
g2=gamma(1+1./beta_vec);
cv=sqrt(g1-g2.^2)./g2;
mu=mean(vec);
cv_es=std(vec)/mean(vec);
[min_difference, array_position] = min(abs(cv - cv_es));
beta=beta_vec(array_position);
eta=(mu/gamma(1/beta+1));
end
代码测试:
使用MATLAB的wblrnd(η,β,10000,1)产生测试数据,使用上述方法进行参数估计。
β | η | β(估计值) | η(估计值) |
---|---|---|---|
4 | 2 | 3.9710 | 2.0077 |
1 | 3 | 1.0040 | 3.0231 |
0.4 | 0.5 | 0.4100 | 0.5440 |