Censored Weibull Distribution 最大似然估计 （结合牛顿法求解）

前言：写这篇博客是因为我前几天偶然读到一篇很有意思的文章， 然后想用自己的实验数据测试一下其分布情况，本以外是一个很简单的工作，但再简单的工作也需要不断的更新学习。在学习过程中，我发现已经有很多博文详述了Weibull分布的应用和意义，然而对于如何进行参数估计的问题，却缺乏更加详细的记录文档。于是无聊而又不务正业的Shecan稍微探索了一下，并把笔记分享给大家。如有错误或者疑问，请给Shecan留言。

Weibull 分布函数 CDF：
$F(x) = 1 - e^{-(\frac{x}{\beta})^{\alpha}}\tag{1}$
这里 $\beta$ 表示 Scale Parameter, $\alpha$ 表示 Shape Parameter。

分布密度 PDF：
$f\left( x\right) =\frac {\alpha } {x}\gamma ^{\alpha }e^{-\gamma ^{\alpha}}\tag{2}$
其中
$\gamma = \frac{x}{\beta}$
censored 似然函数的定义：
$L=\prod _{i=1}^{n}\left( f\left( x_i\right) \right) ^{\delta_i}\left( 1-F\left( x_{i}\right) \right) ^{1-\delta_{i}}\tag{3}$
其中
$\delta _{i}=\begin{cases} 1,\,x\leq threshold\\ 0,\,x>threshold\end{cases}$
把（1），（2）带入（3）两边求对数，并简化求得对数似然函数：
$1og_L=\sum _{i=1}^{n}\left( - \gamma _{i}^{\alpha }+\delta_{i}\alpha ln\gamma _{i}+\delta\ln \frac {\alpha } {x_{i}}\right)\tag{4}$
其中 $\alpha$， $\beta$, 可以通过最大化似然函数求得，对（4）求一阶导数
$\sum_{i=1}^{n}\gamma_i^\alpha-\delta_i=0,\tag{5}$
$\sum _{i=1}^{n}\left( -\gamma _{i}^{\alpha }\ln \gamma _{i}+\delta_i\ln \gamma _{i}+\frac {\delta_i} {\alpha }\right) =0.\tag{6}$
（5），（6）是关于 $\alpha,\,\beta$ 的外生解，无法求出其解析解，因而需要求解以上两个非线性方程组的数值解。求解非线性方程组数值解的方法有很多，比如非线性最小二乘法，Broyden方法等等，但是我们遵循简单即有效的原则，尝试最为简单的Newton-Rapson解的搜索方法。首先，简化（5）得
$\beta =\left( \frac {\sum_{i}^{n}x_{i}^{\alpha }} {\sum_{i}^{n}\delta_{i}}\right) ^{\frac {1} {\alpha }}\tag{7}$
因此，我们只需要通过（6）求解 $\alpha$, 假设
$h(\alpha)=\sum _{i=1}^{n}\left( -\gamma _{i}^{\alpha }\ln \gamma _{i}+\delta_i\ln \gamma _{i}+\frac {\delta_i} {\alpha }\right)\tag{8}$
其中
$\gamma _{i}=\frac {x_{i}} {\beta }=x_i\left( \frac {\sum x_{i}^{\alpha}}{n}\right)^{-\frac {1} {\alpha }}$
and
$\gamma _{i}^{\alpha }=x_{i}^{\alpha}\cdot \frac {n} {\sum x_{i}^\alpha} \tag{9}$
$\ln \gamma _{i}=\ln x_i-\frac {1} {\alpha }\cdot \frac {\sum x_{i}^{\alpha }} {n} \tag{10}$
这里的 $n=\sum \delta_i$ 代表 failed number。

（9）和 （10）代入（8），即得
$h\left( \alpha \right) =n\left( -\frac {\Sigma x ^{\alpha }_iln x_{i}} {\sum x_{i}^{\alpha }}+\frac {1} {\alpha }+\frac {1} {n}\Sigma \delta_i\ln x_{i}\right)\tag{10}$
求导可得
$h'\left( \alpha \right) =n\left( -\frac {1} {\alpha^2}+\frac {\sum x_{i}^{\alpha }\left( ln x_{i}\right) ^{\alpha }} {\sum x_{i}^{\alpha }}-\frac {\left( \sum x_{i}^{\alpha }\ln x_{i}\right) ^{2}} {\left( \sum x_i^{\alpha }\right) ^{2}}\right)$
因而，我们通过牛顿法迭代搜索非线性方程 $h(\alpha)=0$ 的解。
$\alpha_{k+1}=\alpha_{k}-\frac{h\left(\alpha_{k}\right)}{h^{\prime}\left(\alpha_{k}\right)}\tag{11}$

求解出 $\alpha$ 后，通过 （7）即可求解 $\beta$ 。

Weibull MLE 程序实现：

推导过程写完了，怎么用程序实现呢？从多年来失败的研究经验Shecan总结了一个道理：不要重复造轮子。（其实学习的时候重复造下轮子也有好处。） Shecan 发现Python有一个比较有意思的Package。 叫做Weibull，基本可以满足你大部分的研究需求。具体怎么使用，请参考下面的参考文献。但是，有时候我们依旧需要对模型进行改进，比如说两参数的Weibull分布改成三参数的Weibull分布。为满足这种需求，Shecan找出了Weibull Package MLE部分的源程序，并改写成了面向过程的编程方式，具体如下：

你可以把下面的代码复制到 Jupyter Notebook：

注明：程序来自于Weibull Package:

[https://weibull.readthedocs.io/en/latest/examples.html](Weibull Package Doc)

import numpy as np
import pandas as pd

测试数据：

fail_times = [ 9402.7, 6082.4, 13367.2, 10644.6, 8632.0, 3043.4, 1034.5, 2550.9, 2550.9, 3637.1]
suspended = [True, True, True, True, True,False, False, False, True, True]
data = pd.DataFrame({'data': fail_times, 'susp':suspended})

MLE Calibration: 通过以上的推导过程可以很容易的读懂这段代码。

# filter the failed samples and extract values
df_failed = data[data.susp == False].copy() 
dtf_failed = df_failed["data"].values
df_failed["ln_x_div_r"] = df_failed.apply(lambda s: np.log(s['data'])/len(df_failed), axis=1)
# extract data of all the values
dtf_all = data['data'].values
# use Newton-Rhapson method for estimating the shape parameter

# give initial value for the shape paramter:
shape = (((6.0 / np.pi ** 2)
            * (np.sum(np.log(dtf_all) ** 2)
            - ((np.sum(np.log(dtf_all))) ** 2) / dtf_all.size))
            / (dtf_all.size - 1)) ** -0.5

# 10 iterations of the newton-rhapson method
for i in range(1, 11):
    a = np.sum(np.log(dtf_failed) * 1.0) / dtf_failed.size
    b = np.sum(dtf_all ** shape)
    c = np.sum((dtf_all ** shape) * np.log(dtf_all))
    h = np.sum((dtf_all ** shape) * (np.log(dtf_all)) ** 2)

    shape = shape + (a + (1.0 / shape) - (c / b)) / ((1.0 / shape ** 2) + ((b * h) - c ** 2) / b ** 2)

shape = max(shape, 0.005)
scale = (np.sum((dtf_all ** shape) / len(df_failed))) ** (1 / shape)

# Print the results
print(shape,scale)

题外话：小时候老师教育我们如果要成为一个成功的人，头悬梁，锥刺骨，卧薪尝胆，这些典故都很励志。可是长大后发现了一个残酷的真相，成功人士万里挑一，而千千万万的普通人虽然没有那么成功，却也在以一己之力改变世界。比如说战疫前线的医生和护士们，比如说用自己的才能为这次疫情做贡献的各路神仙。经历了这些事，这些人，Shecan开始反思，与其教育自己的孩子成为一个成功的人，不如教育他们尽一己之力做一个对社会有贡献的人，或许这样的人生会更自信，更快乐，并获得更高的自我效能感。这也算是Shecan开始做笔记，并分享笔记的原因吧。