用 MLE 和 Newton-Raphson 估计 Weibull

Posted 2023-02-19

【中文标题】用 MLE 和 Newton-Raphson 估计 Weibull

【英文标题】：Estimating Weibull with MLE and Newton-Raphson

【发布时间】：2021-06-20 15:54:54

【问题描述】：

我一直在尝试使用牛顿法估计二参数 Weibull 分布。在阅读有关使用 Newton-Raphson 算法的一些信息时，我发现理解某些方面具有挑战性。

我尝试在 Python 中实现它，并且我认为我的方法没有错。但是由于我一直在努力理解算法本身，所以我认为我遗漏了一些东西。我的代码运行，问题是它没有找到正确的估计（1.9 和 13.6）：

#data input in  the Weibull dist.
t = np.array(list(range(1, 10)))
t = np.delete(t,[0])

#calculating the first and second partial derivative of Weibull log-likelihood function
def gradient(a,b): 
    for i in t: 
        grad_a = np.array(-10*b/a + b/a*np.sum((i/a)**b),dtype = np.float)
        grad_b = np.array(10/b - 10*(math.log(a)) + np.sum(math.log(i)) - np.sum(((i/a)**b)*math.log(i/a)),np.float)     
        grad_matrix = np.array([grad_a, grad_b])
    return grad_matrix
    
def hessian(a,b): 
    for i in t: 
        hess_a = np.array((10*b/a**2 + (b*(b+1)/a**2)*np.sum((i/a)**b)),np.float)
        hess_b = np.array(10/b**2 + np.sum(((i/a)**b) * (math.log(i/a))**2),np.float)
        hessians = np.array([hess_a, hess_b]) 
    return hessians  

#Newton-Raphson
iters = 0     
a0, b0 = 5,15

while iters < 350:  
    if hessian(a0,b0).any() == 0.0:
        print('Divide by zero error!') 
    else:
        a = a0 - gradient(a0,b0)[0]/hessian(a0,b0)[0]
        b = b0 - gradient(a0,b0)[1]/hessian(a0,b0)[1]    
        print('Iteration-%d, a = %0.6f, b= %0.6f, e1 = %0.6f, e2 = %0.6f' % (iters, a,b,a-a0,b-b0))    
    if math.fabs(a-a0) >0.001 or math.fabs(b-b0) >0.001:
        a0,b0 = a,b
        iters = iters +1
    else: 
        break
print(a,b)
print(iters)    

**Output:**             
Iteration-0, a = 4.687992, b= 16.706941, e1 = -0.312008, e2 = 1.706941          
Iteration-1, a = 4.423289, b= 18.240714, e1 = -0.264703, e2 = 1.533773                
Iteration-2, a = 4.193403, b= 19.648545, e1 = -0.229886, e2 = 1.407831     

     

以此类推，每次迭代都离第二个参数 (b) 的正确估计越来越远。

威布尔 pdf: http://www.iosrjournals.org/iosr-jm/papers/Vol12-issue6/Version-1/E1206013842.pdf

【问题讨论】：

你能给出你的 2-param Weibull 分布方程吗？我想检查你的梯度和粗麻布。顺便说一句，在我看来，您只是在 for 循环中覆盖了 grad_a 和 grad_b，而不是使用 +=。但是，如果没有确切的符号，我无法轻松验证您的代码。牛顿部分似乎还可以。

@flow_me_over，非常感谢您确认 NR 至少看起来还不错！我使用了以下 Weibull pdf：f(t; a, b) = b/a * (t/a)^(b-1)*exp-(t/a)^b。它对应于等式。 （3.1）在我编辑的帖子中附加的论文中，我也从中获取了渐变和粗麻布。导数取自 Weibull pdf 的对数似然。

@flow_me_over，问题是我使用连续 Weibull pdf 导出导数，而我的 t 是离散的...

【参考方案1】：

好的。所以首先我提一下，你使用的论文不清楚，令我惊讶的是，这项工作能够进入期刊。其次，您声明您的输入数据“t”，即论文中的“x”，是从 0 到 9 的数字列表？我在论文中找不到这个，但我假设这是正确的。

下面我更新了您的渐变函数，该函数非常冗长且难以阅读。我已经使用 numpy 为你矢量化了它。看看你是否明白。

您的粗麻布不正确。我相信论文中的二阶导数中有一些错误的迹象，因此在你的。也许再过一遍它们的推导？尽管如此，无论符号如何变化，您的 Hessian 都没有明确定义。一个 2x2 Hessian 矩阵包含对角线上的二阶导数 d^2 logL / da^2 和 d^2 logL /db^2，以及非对角线上的导数 d^2 log L /da db（位置 (1,2 ) 和 (2,1) 在矩阵中)。我已经调整了你的代码，但并不是我没有纠正可能错误的迹象。

总之，您可能需要根据 Hessian 更改再次检查您的 NR 代码，并创建一个 while 循环，以确保算法在您达到容差水平后停止。

#data input in the Weibull dist.
t = np.arange(0,10) # goes from 0 to 9, so N=10
N = len(t)

#calculating the first and second partial derivative of Weibull log-likelihood 
#function
def gradient(a,b): 
    grad_a = -N*b/a + b/a*np.sum(np.power(t/a,b))
    grad_b = N/b - N*np.log(a) + np.sum(np.log(t)) - np.sum(np.power(t/a,b) * np.log(t/a)))      

    return np.array([grad_a, grad_b])

def hessian(a,b):  
    hess_aa = N*b/a**2 + (b*(b+1)/a**2)*np.sum(np.power(t/a,b))
    hess_bb = N/b**2 + np.sum(np.power(t/a,b) * np.power(np.log(t/a),2))
    hess_ab = ....
    hess_ba = hess_ab
    hessians = np.array([[hess_aa, hess_ab],[hess_ba, hess_bb]]) 
    
    return hessians  

希望这些 cmets 能进一步帮助您！请注意，Python 有一些库可以在数值上找到对数似然的最佳值。参见例如 scipy 库，函数“最小化”。

【讨论】：

非常感谢！是的，确实，我手动计算了粗麻布矩阵的元素，现在我的 NR 开始了！不重新检查论文中的衍生物是我的天真，所以它发表了很奇怪///