# Newton-Raphson for estimating Weibull distribution does not converge

I've been trying to estimate the two-parameter (a,b) Weibull distribution (loc. param. = 0).

$$f(t;a,b)=\frac{b}{a}\left(\frac{t}{a}\right)^{b-1}\exp(- \left(\frac{t}{a}\right)^b)$$

To find the estimates I am using maximum likelihood and NR method.

t = np.array(list(range(1, 1000))) #data input
grad_b = 10/b - 10*(np.log(a)) + np.sum(np.log(t)) - np.sum(((t/a)**b)*np.log(t/a))

def hessian(a,b):
hess_a = 10*b/a**2 + (b*(b+1)/a**2)*np.sum((t/a)**b)
hess_b = 10/(b**2) + np.sum((t/a)**b* (np.log(t/a))**2)
hessians = np.array([hess_a, hess_b])
return hessians

iters = 0
a0, b0 = 2,6 #initial condition

while iters < 4:
if hessian(a0,b0).any() == 0.0:
print('Divide by zero error!')
else:
print('Iteration-%d, a = %0.6f, b= %0.6f, e1 = %0.6f, e2 = %0.6f' % (iters, a,b,a-a0,b-b0))
break
else:
a0,b0 = a,b
iters = iters +1
print(a,b)
print(iters)


Output:
Iteration-0, a = 1.714391, b= 6.656262, e1 = -0.285609, e2 = 0.656262
Iteration-1, a = 1.490482, b= 7.249362, e1 = -0.223909, e2 = 0.593099
Iteration-2, a = 1.309805, b= 7.794966, e1 = -0.180677, e2 = 0.545604
Iteration-3, a = 1.160878, b= 8.303167, e1 = -0.148927, e2 = 0.508201

a goes to zero, b to infinity, neither of which is a correct estimation. It feels like I am doing something wrong in defining gradients and hessians, though these seem to be correct first and second derivatives of the log-likelihood function.

• I'm curious what the factors of 10 in your gradient and Hessian are for. Mar 24 at 23:56
• @Aruralreader, this is just the size of t (number of data points), because I am taking gradient and hessian from the log-likelihood function Mar 25 at 0:07