I'm currently trying to estimate parameters of a distribution with the mle method in Python. These are the derivative of my loglikelihood function: Loglikelihood Partial Derivatives
As you can see it has a quite complicated formula and I would need to derive the roots numerically (I want to use the Newton Raphson method) this is how I have written my code:
from math import cos, sin, pi, exp, sqrt, log
from scipy.special import psi
import autograd as ad
from autograd import grad, jacobian
import autograd.numpy as np
import sympy as sy
x = [0.1, 0.2, 1, 1, 1, 1, 1, 2, 3, 6, 7, 11, 12, 18, 18, 18, 18, 18, \
21, 32, 36, 40, 45, 46, 47, 50, 55, 60, 63, 63, 67, 67, 67, 67, 72, 75, 79, \
82, 82, 83, 84, 84, 84, 85, 85, 85, 85, 85, 86, 86]
p = [2.5, 0.00003, 0.1, 0.02, 0.1]
p0 = p[0]
p1 = p[1]
p2 = p[2]
p3 = p[3]
p4 = p[4]
def Fs():
n = len(x)
sum1 = 0
sum2 = 0
sum3 = 0
sum4 = 0
sum5 = 0
sum6 = 0
sum7 = 0
sum8 = 0
sum9 = 0
sum10 = 0
for i in range(1,n):
z = exp(-(p1 / p2) * (exp(p2 * x[i]) - 1))
a = 1 - z
b = a ** p0
c = exp(p2 * x[i])
sum1 = sum1 + log(a)
sum2 = sum2 + b * log(a) / log(1 - b)
sum3 = sum3 + c
sum4 = sum4 + (c - 1) * z / a
sum5 = sum5 + (c - 1) * z * (a ** (p0 - 1)) / (1 - b)
sum6 = sum6 + x[i]
sum7 = sum7 + c * (1 - p2 * x[i])
sum8 = sum8 + (p2 * x[i] * c - c + 1) * z / a
sum9 = sum9 + (p2 * x[i] * c - c + 1) * z * \
(a ** (p0 - 1)) / (1 - b)
sum10 = sum10 + log(1 - b)
f1 = n / p0 + p3 * sum1 - (p4 - 1) * sum2
f2 = n / p1 + n / p2 - sum3 / p2 + sum4 * (p0 * p3 - 1) / \
p2 + sum5 * p0 * (p4 - 1) / p2
f3 = -n * p1 / (p2 ** 2) + sum6 + p1 * sum7 / (p2 ** 2) - \
p1 * (p3 * p0 - 1) * sum8 / (p2 ** 2) + \
p0 * p1 * (p4 - 1) * sum9 / (p2 ** 2)
f4 = -n * (psi(p3) + psi(p3 + p4)) + p0 * sum1
f5 = -n * (psi(p4) + psi(p3 + p4)) + p0 * sum10
return np.matrix([[f1], [f2], [f3], [f4], [f5]])
So the problem here is, I would like to get the Jacobian matrix for the NR method using jacobian(). However I realize that if I run my Fs function the output will be in numbers not in the form of a polynom, therefore making the Jacobian function unable to perform. I was wondering can anybody help me to have my functions as polynom instead of values? Are there any way I could represent the sum in my functions without going in a loop? Or any suggestions how I might tackle this problem? Thank you in advance.
