韦伯分布(Weibull)参数矩估计MATLAB实现

二参数韦伯分布概率密度函数

$$f(x)=\frac{\beta }{\eta }{\left(\frac{x}{\eta }\right)}^{\beta -1}{e}^{-{\left(\frac{x}{\eta }\right)}^{\beta }},\beta >0,\eta >0,x\ge \gamma \ge 0$$

其中,$\beta$是形状参数，$\eta$是尺度参数。

$\beta 和\eta$的矩估计

1. Weibull分布的k阶矩：

$${\mu }_k={\left(\frac{1}{{\eta }^{\beta }}\right)}^{-\frac{k}{\beta }}\Gamma \left(1+\frac{k}{\beta }\right)$$

其中，$\Gamma$是gamma函数：
$$\Gamma (s)={\int }_0^{\infty }{x}^{s-1}{e}^{-x}dx,(s>0)$$
推导过程：
$$\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 }}dx\\ ={\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}dx\\ ={\eta }^k\Gamma \left(\frac{k}{\beta }+1\right)\end{array}$$

2. 因此可以计算出一阶矩和二阶矩

$$\begin{array}{rl}{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{array}$$

让m2除以m1的平方，我们可以得到：
$$\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)}$$
这正好是变异系数的平方，上式开方就是变异系数：
$$CV=\frac{\sqrt{\Gamma \left(1+\frac{2}{\beta }\right)-{\Gamma }^2\left(1+\frac{1}{\beta }\right)}}{\Gamma \left(1+\frac{1}{\beta }\right)}$$

3. 我们可以借助变异系数CV来估计$\beta$

1. 首先我们定义一组候选$\beta$，如 { β ∣ 0.1 < β < 10 } \{\beta|0.1<\beta<10\} {β∣0.1<β<10}的数组，根据这一系列候选$\beta$计算出对应的$CV$数组
2. 然后根据数据的均值和方差计算数据的变异系数$\widehat{CV}$
3. 找出$CV$数组中和$\widehat{CV}$最接近的值，其对应的$\beta$就是我们需要的估计值

4. 尺度参数$\eta$可以利用下面的公式估计

η ^ = { x ˉ / Γ [ ( 1 / β ^ ) + 1 ] } \hat{\eta}=\{\bar{x} / \Gamma[(1 / \hat{\beta})+1]\} η^​={xˉ/Γ[(1/β^​)+1]}

其中$\bar{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

参考文献

AlFawzan, Ma. Methods for estimating the parameters of the Weibull distribution[J]. Fawzan, 2000(10).