用最大似然估计和牛顿-拉夫森估计威布尔

Question

我一直在尝试用牛顿法估计两参数威布尔分布。当我读到一些关于使用牛顿-拉夫森算法的文章时，我发现理解某些方面是很有挑战性的。

我试着用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对它进行了矢量化。看看你是不是明白了。

你的黑森语不正确。我相信在论文中的二阶导数中有一些错误的信号，因此在你的论文中也是如此。也许再看一遍它们的推导过程？尽管如此，不管符号如何变化，你的黑森语都不是很好的定义。2x2Hessian矩阵包含对角线上的二阶导数d^2 logL / da ^2和d^2 logL /db^2，以及非对角线上的导数d^2log 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提供了一些库，可以用数值方法找到对数似然率的最优值。例如，参见scipy库，函数'minimize‘。