使用 MLE 和 Newton-Raphson 估计 Weibull

Question

我一直在尝试用牛顿法估计二参数威布尔分布。当我阅读一些有关使用 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

Answer 1

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

下面我更新了你的渐变函数，它非常冗长且难以阅读。我已经使用 numpy 为您对其进行了矢量化。看懂了没。

您的 Hessian 矩阵不正确。我相信论文中的二阶导数有一些错误的迹象，因此在你的。也许再回顾一下它们的推导？尽管如此，无论符号如何变化，您的 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

希望这些评论对您有进一步的帮助！请注意，Python 有一些库可以在数值上找到对数似然的最优值。例如参见 [=25=] 库，函数 'minimize'.

使用 MLE 和 Newton-Raphson 估计 Weibull

Estimating Weibull with MLE and Newton-Raphson

python

derivative

newtons-method

mle

weibull