三参数Weibull分布 $W(\beta ,\eta ,\gamma )$的概率密度分布函数为
$$f(t;\eta ,\beta ,\gamma )=\frac{\beta }{\eta }(\frac{t-\gamma }{\eta }{)}^{\beta -1}{e}^{-(\frac{t-\gamma }{\eta }{)}^{\beta }}\cdots \cdots （1）$$
式中， $\eta$、 $\beta$、和 $\gamma$分别称为尺度参数、形状参数和位置参数。当 $\gamma =0$时，模型退化为两参数威布尔分布。

  三参数Weibull分布$W(\beta ,\eta ,\gamma )$的累积分布函数为
$$F(t)=1-e^{-(\frac{t-\gamma }{\eta }{)}^{\beta }}\cdots \cdots （2）$$

  式中，形状参数$\beta$具有域$(0,\infty )$，它被称为$Weibull-$斜率，因为它给出了在$Weibull-$概率纸上的斜率。

一、相关系数法

  相关系数法（Correlation coefficient method，CCM）是当样本的相关系数方程式（4）取得最大值时的位置参数即为估计的位置参数。再使用最小二乘法式（7）求得形状参数和尺度参数。
将$F(t)$取两次对数得到
$$ln(-ln(1-F))=\beta ln(t-\gamma )-\beta ln\eta \cdots \cdots （3）$$
  则$ln(-ln(1-F))$与$ln(t-\gamma )$呈线性关系，形如$y=\beta x+b$。这里，$y=ln(-ln(1-F))$，$x=ln(t-\gamma )$，$b=-\beta ln\eta$。样本的相关系数为
$$R(\gamma )=\frac{{\sum }_{i=1}^n{x}_i{y}_i-n\bar{x}\cdot \bar{y}}{\sqrt{({\sum }_{i=1}^n{x}_i^2-n{\bar{x}}^2)({\sum }_{i=1}^n{y}_i^2-n{\bar{y}}^2)}}\cdots \cdots （4）$$
式中，$x_i=ln(t_i-\gamma )$，$\bar{x}={\sum }_{i=1}^n{x}_i/n$，$y_i=ln(-ln(1-F_i))$，$\bar{y}={\sum }_{i=1}^n{y}_i/n$。式（3）中$F$为累积分布函数$F(t)$，可以使用中位秩来估计
$$\hat{F}(t_i)\approx \frac{i-0.3}{n+0.4},(i=1,2,\cdots ,n)\cdots \cdots （5）$$
式中，$i$是按升序排列的寿命样本的序号；$n$是样本量。
  在实际应用中位置参数添加一个条件$0\le \gamma <t_{(1)}$是十分合理的。然后，位置参数$\gamma$的估计由式（6）给出
$${\hat{\gamma }}_{CCM}=argmaxR(\gamma )\cdots \cdots （6）$$
  获得位置参数估计后，将原来的样本减去位置参数后得到新的样本使用最小二乘方程式（7）可得到形状参数和尺度参数
$$\{\begin{array}{c}{\hat{\beta }}_{CCM}=\frac{{\sum }_{i=1}^n{x}_i{y}_i-n\bar{x}\cdot \bar{y}}{{\sum }_{i=1}^n{x}_i^2-n{\bar{x}}^2}\\ {\hat{\eta }}_{CCM}=exp(-\frac{\bar{y}-\hat{\beta }\bar{x}}{\hat{\beta }})\end{array}\cdots \cdots （7）$$

二、最小二乘法

  最小二乘法（Least squares method，LSM）是将累积分布函数的偏差平方和作为性能指标，累积分布函数的偏差平方和最小时，则$\beta$，$\eta$，$\gamma$为求得的估计值
$$({\hat{\beta }}_{LSM},{\hat{\eta }}_{LSM},{\hat{\gamma }}_{LSM})=arg\underset{(\beta ,\eta ,\gamma )}{\min }\sum _{i=1}^n[F(t_i)-F(t){]}^2\cdots \cdots （8）$$
其中，累积分布函数可以用中位秩式（5）估计。

三、最小偏差法

  最小偏差法（Minimum discrepancy method，MDM）由谢里阳等提出。通过转换累积分布函数，构建了一个从随机变量、其对应的累积分布概率、形状参数和位置参数到尺度参数的映射。当形状参数和位置参数出现误差时，这种映射产生的尺度参数取决于随机变量和累积分布概率。给定一组样本值，形状/位置参数的较大误差会使得估计出的尺度参数值之间产生较大的差异。通过将尺度参数之间的差异最小化，可估计出正确的形状参数和位置参数。
  由式（2）可转换为
$$\eta =\frac{t-\gamma }{(-ln(1-F){)}^{1/\beta }}\cdots \cdots （9）$$
式中，$t$为寿命数据，其中累积分布函数可以用中位秩式（5）估计。若形状参数准确，则每一个寿命数据与其对应的失效概率计算出的尺度参数应该相等，所以当计算的多个尺度参数的标准差最小时，此时的形状参数$\beta$、位置参数$\gamma$为估计的形状参数${\hat{\beta }}_{MDM}$和估计的位置参数${\hat{\gamma }}_{MDM}$，则形状参数和位置参数的计算公式为
$$({\hat{\beta }}_{MDM},{\hat{\gamma }}_{MDM})=arg\underset{(\beta ,\gamma )}{\min }(std(\frac{t_i-\gamma }{(-ln(1-F_i){)}^{1/\beta }}))\cdots \cdots （10）$$
  为了消除样本值对尺度参数估计的不确定性影响，可将平均值用于威布尔尺度参数估计
$${\hat{\eta }}_{MDM}=\frac{1}{n}\sum _{i=1}^n\frac{t_i-{\gamma }_{MDM}}{(-ln(1-F_i){)}^{1/{\hat{\beta }}_{MDM}}}\cdots \cdots （11）$$

四、极大似然估计法

• 这是一种广泛应用的参数估计方法，在Weibull分布中也同样适用。
• 通过构建似然函数，即将观测数据代入Weibull分布的概率密度函数，然后最大化这个似然函数来求解形状参数和尺度参数。
• 对数似然函数是似然函数的自然对数，通过对其求偏导并令导数等于零，可以求解出参数的估计值。

  对于容量为$n$的完全样本数据$t_1\le t_2\le \cdots \le t_n$，Weibull分布的对数似然函数为
$$lnL(t_1,\cdots ,t_n;\eta ,\beta ,\gamma )=\sum _{i=1}^{n}ln(\frac{\beta }{{\eta }^{\beta }}(t_i-\gamma {)}^{\beta -1}exp(-(x_i-\gamma {)}^{\beta }/{\eta }^{\beta }))\cdots \cdots （12）$$
  对数似然函数对各参数求偏导数 ，得方程组
$$\frac{\partial lnL}{\partial \eta }=-\frac{n}{\eta }-\frac{n(\beta -1)}{\eta }+\frac{\beta }{\eta }\sum _{i=1}^n(\frac{t_i-\gamma }{\eta }{)}^{\beta }=0\cdots \cdots （13）$$
$$\frac{\partial lnL}{\partial \beta }=\frac{n}{\beta }+\sum _{i=1}^{n}ln\frac{t_i-\gamma }{\eta }-\sum _{i=1}^n(\frac{t_i-\gamma }{\eta }{)}^{\beta }ln(\frac{t_i-\gamma }{\eta })=0\cdots \cdots （14）$$
$$\frac{\partial lnL}{\partial \gamma }=(1-\beta )\sum _{i=1}^n\frac{1}{t_i-\gamma }+\frac{\beta }{\eta }\sum _{i=1}^n(\frac{t_i-\gamma }{\eta }{)}^{\beta -1}=0\cdots \cdots （15）$$
  上述方程组即为求解威布尔分布参数估计的似然方程组，解方程组（13）～（15）即可获得 3 个参数的估计值。但 由于方程组（13）～（15）无代数解 ，求解十分复杂 ，一般是通过计算机语言如 c 语言、M atlab 等编程求解 。
  以下通过直接寻求对数似然 函数的极大值求解各参数。
  由$\partial L(t_1,...,t_n;\eta ,\beta ,\gamma )/\partial \eta =0$可求得
$${\eta }^{\beta }=\frac{1}{n}\sum _{i=1}^n(t_i-\gamma {)}^{\beta }\cdots \cdots （16）$$
把（16）式代入（12）式 ，从（12）式中消去尺度参数$\eta$可得
$$lnL(t_1,\cdots ,t_n;\beta ,\gamma )=nln\beta -nln\left[\frac{1}{n}\sum _{i=1}^n(t_i-\gamma {)}^{\beta }\right]+(\beta -1)\sum _{i=1}^{n}ln(t_i-\gamma )-n\cdots \cdots （17）$$
  根据极大似然原理，使（17）式取得极大值的$\hat{\beta }$，$\hat{\gamma }$即为所求的形状参数、位置参数的估计值，然后再根据（16）式即可估计出尺度参数$\hat{\eta }$。由三参数Weibull分布的物理意义可知，一般有$0<\beta <10$，$0\le \gamma <x_1$，这样就把求解对数似然方程组的问题变为了求解有约束条件的极值问题，使问题得到了大大简化，而且（17）式的极值可在MS EXCEL上使用“规划求解”功能直接求解，方便了一般工程技术人员使用。极大似然估计法是一种常用的参数估计方法，精度较高，且适用于包括有中途撤出试验的各种截尾试验场合。

五、矩法估计

  若设$g_i(\beta )=\Gamma (1+\frac{i}{\beta }),i=1,2,3$，则三参数Weibull分布的数学期望$E(X)$、标准差$\sigma (X)$、偏度$B(X)$及其估计值样本均值$\bar{x}$，样本标准差$S$，样本偏度$B_0$分别是
$$E(X)=\eta \cdot g_1(\beta )+\gamma \cdots \cdots （18）$$
$$\sigma (X)=\eta \sqrt{g_2(\beta )-g_1^2(\beta )}\cdots \cdots （19）$$
$$B(X)=\frac{g_3(\beta )-3g_2(\beta )g_1(\beta )+2g_1^3(\beta )}{{\left[g_2(\beta )-g_1^2(\beta )\right]}^{3/2}}\cdots \cdots （20）$$
$$\bar{x}=\frac{1}{n}\sum _{i=1}^n{x}_i\cdots \cdots （21）$$
$$S=\sqrt{\frac{1}{n-1}\sum _{i=1}^n(x_i-\bar{x}{)}^2}\cdots \cdots （22）$$
$$B_0=\frac{n}{(n-1)(n-2)S^3}\sum _{i=1}^n(x_i-\bar{x}{)}^3\cdots \cdots （23）$$
  利用（23）式求出$B_0$后作为$B(X)$的估计值，代入（20）式可求出β的估计值$\hat{\beta }$，然后用（19）式和（22）式获得$\eta$的估计值$\hat{\eta }$，再用（18）式和（21）式获得$\gamma$的估计值$\hat{\gamma }$。
  矩法估计的基本思想是用试验样本的各阶矩估计母体的各阶矩，并据此估计其他参数。矩法估计算法简单，有专用数表可查，使用方便。该法在小样本时精度不高，且仅适用于完全样本的场合。

六、概率权重矩法

  概率权重矩法不需要迭代计算 ，但当样本容量较小时，该方法精度较差。
  试验样本概率权重矩的估计值为
$$M_{1,0,k}=\frac{1}{n}\sum _{i=1}^n{x}_i(1-\frac{i-0.35}{n}{)}^k,k=0,1,3\cdots \cdots （24）$$
用概率权重矩表示的Weibull分布参数的估计值为
$$\hat{\eta }=\frac{M_{1,0,0}-\gamma }{\Gamma \left[ln(\frac{M_{1,0,0}-2M_{1,0,1}}{M_{1,0,1}-2M_{1,0,3}})/ln2\right]}\cdots \cdots （25）$$
$$\hat{\beta }=\frac{ln2}{ln(\frac{M_{1,0,0}-2M_{1,0,1}}{2M_{1,0,1}-4M_{1,0,3}})}\cdots \cdots （26）$$
$$\hat{\gamma }=\frac{4(M_{1,0,3}{M}_{1,0,0}-M_{1,0,1}^2)}{4M_{1,0,3}+M_{1,0,0}-4M_{1,0,1}}\cdots \cdots （27）$$
  概率权重矩法的精度与样本概率权重矩的计算方法有很大的关系，同时与数据的分散性有较大的关系，有时无法获得满足要求的解。在各种估计方法中，该法估计出的位置参数较其他方法小。

七、双线性回归法

  变换（2）式，并令$t_0={\eta }^{\beta }$，可得到两个方程
$$ln\left[-ln(1-F(t))\right]=\beta ln(t-\gamma )-ln{t}_0\cdots \cdots （28）$$
$${\left[-ln(1-F(t))\right]}^{1/\beta }=\frac{1}{\eta }t-\frac{\gamma }{\eta }\cdots \cdots （29）$$
  可以证明上述两式线性无关，每个方程单独进行最小二乘估计，合并整理后得到
$$\beta =\frac{ln(t-\gamma )\cdot ln\left[-ln(1-F(t))\right]-ln(t-\gamma )\cdot ln\left[-ln(1-F(t))\right]}{l{n}^2(t-\gamma )-ln(t-\gamma {)}^2}\cdots \cdots （30）$$
$$\gamma =\bar{t}-\frac{{\left[-ln(1-F(t))\right]}^{1/\beta }\cdot (t^2-{\bar{t}}^2)}{t\cdot {\left[-ln(1-F(t))\right]}^{1/\beta }-\bar{t}\cdot {\left[-ln(1-F(t))\right]}^{1/\beta }}\cdots \cdots （31）$$
$${\eta }^{\beta }=exp\left[\beta \cdot ln(t-\gamma )-ln(-ln(1-F(t)))\right]\cdots \cdots （32）$$
  任给$\beta$的一个初始值，由式得（31）式得$\gamma$的一个初值，代入（30）式得$\beta$的一个估计值$\hat{\beta }$值，若$\beta$和$\hat{\beta }$差值足够小，则$\hat{\beta }$即为所求的形状参数，否则用$2\beta -\hat{\beta }$代替$\beta$，继续用（31）式计算$\gamma$，再代入（30）式计算$\hat{\beta }$，如此反复直至满足精度要求。最后用（31）式计算尺度参数$\hat{\eta }$。
  双线性回归法是一种精度较高的方法，但在迭代过程中，有些样本会出现负数取对数的现象，使得在MS EXCEL上无法使用宏功能求解，这种情况下可使用其他方法代替。

八、灰色模型法

  灰色模型法基于邓聚龙提出的灰色系统原理是一种较新的参数估计方法。该法在数据量较小时就可获得较高的估计精度。灰色模型法也无需迭代计算，用最小二乘法即可获得3个参数的估计值，在各种样本容量下均可获得较好的精度。
   不同容量的样本实例计算表明，随着样本容量的增加，各估计方法之间的差异越来越小，在样本容量较大时，各种估计方法均可使用；在样本容量较小时，各种估计方法的差异较大，灰色模型法具有较高的精度。
  （2）式也可表示为
$$t=\gamma +\eta \cdot exp(\frac{ln[-ln(1-F(x))]}{\beta })\cdots \cdots （33）$$
  令$k_i=ln[-ln(1-F(t_i))],i=1,2,\cdots ,n$，并记$\eta =c$，$1/m=-a$，$\gamma =b$，视$(t_i,k_i)$ 为一时间序列，则（33）式可转化为
$$t_i=c\cdot exp(-a{k}_i)+b\cdots \cdots （34）$$
  灰色模型$GM(1,1)$的微分方程为
$$\frac{d\hat{t}(k)}{dk}+a\hat{t}(k)=u(k\in R)\cdots \cdots （35）$$
  灰色模型$GM(1,1)$的时间响应模型为
$$\hat{t}(k)=c\cdot exp(-ak)+\frac{u}{a}\cdots \cdots （36）$$
  （34）式和（36）式具有相同的形式，因此可用灰色模型对参数$a$，$u$，$c$进行估计，进而得到$\beta$，$\gamma$，$\eta$的估计值。灰色模型的参数为
$${\left[\begin{array}{cc}a& u\end{array}\right]}^T=(B^{T}B{)}^{-1}{B}^T{Y}_N\cdots \cdots （37）$$
式中，$B=\left[\begin{array}{cc}-(t_1+t_2)/2& 1\\ ⋮& ⋮\\ -(t_{n-1}+t_n)/2& 1\end{array}\right]$，$Y_N={\left[\begin{array}{ccc}\frac{t_2-t_1}{k_2-k_1}& \cdots & \frac{t_n-t_{n-1}}{k_n-k_{n-1}}\end{array}\right]}^T$。
$${\left[\begin{array}{cc}b& c\end{array}\right]}^T=(D^{T}D{)}^{-1}{D}^{T}X\cdots \cdots （38）$$
式中，$D={\left[\begin{array}{ccc}1& \cdots & 1\\ exp(-a{k}_1)& \cdots & exp(-a{k}_n)\end{array}\right]}^{{}^T}$，$X={\left[\begin{array}{ccc}{x}_1& \cdots & x_n\end{array}\right]}^T$。
  由（37）式得到$a$和$u$的估计值，对比（34）式和（36）式即知$\hat{m}=-a^{-1}$，$\hat{\gamma }=b=u/a$。由（38）式得到$c$即得$\hat{\eta }=c$，（38）式算出的$b$可以作为$y$的一个优化值。