I want to run simple Maximum Likelihood estimation in python. I want to try it by using Scipy.optimize.minimize
in python. first I'll explain my model so you can figure out what is going to happen.
by MLE I want to estimate best value for 2 variables that maximizes my objective function!
and what are these variables and what is objective function? in picture bellow you can see my objective function! and my variables are 'Betaa' and 'Kasi'.
Now you find out what is my objective function and what are my variables(except item
in objective function. I'll explain it! please keep reading). I want to implement MLE on dataset(click to download it from my google drive) that you can see in picture bellow:
it's pretty simple dataset! contains 501 rows and 1 column. now that you saw my data, it's time for explaining how MLE is gonna work on it!
for each value in 'Loss' column, if the value is >160, i should calculate Target value by objective function(that i talked about it in first section) and if the value is <160, Target for that will be 0! see picture bellow then you will find out what I'm precisely talking about(also you'll find out what is item
that i talked about it)!
MLE is gonna maximize summation of values in Target
by finding best values for 'Betaa' and 'Kasi'!
so here is script that have been written for this purpose:
from scipy.optimize import minimize
import pandas as pd
import numpy as np
My_DataFrame = pd.read_excel("<FILE_PATH_IN_YOUR_MACHINE>\\Losses.xlsx")
# i did put link of "Losses.xlsx" file in the model explaination section
# so you can download it from my google drive.
loss = My_DataFrame["Loss"]
def obj_function(x):
k = x[0]
b = x[1]
target = np.array([])
for iloss in loss:
if iloss>160:
t = np.log((1/b)*(1+(k*(iloss-160)/b)**((-1/k)-1)))
target = np.append(target , t)
# multiplied by (-1) so now we have maximization
return -np.sum(target)
MLE_Model = minimize(obj_function , [0.7,20] )
print(MLE_Model)
and results of this code:
t = np.log((1/b)(1+(k(iloss-160)/b)**((-1/k)-1)))
fun: nan
hess_inv: array([[1, 0], [0, 1]])
jac: array([nan, nan])
message: 'Desired error not necessarily achieved due to precision loss.'
nfev: 126
nit: 1
njev: 42
status: 2
success:False
x: array([-1033.52630224, 25.32290945])
as you can see in results, success: False
and optimal solution didn't found!
and now i don't know how to fix this issue and how to solve this MLE problem.
any help would be appreciated.
As already mentioned in the comments, there are multiple runtime warnings due to overflow when evaluating the power. Looking closely, you can observe that your objective function is wrong. According to your picture, it should be
(1 + ( k*(iloss-160)/b) )**( (-1/k) - 1 )
instead of
(1 + ( k*(iloss-160)/b )**( (-1/k) - 1 )
Furthermore, you can avoid the calculation of the power and rewrite your objective by using the logarithm rules:
np.log(1/b) + ((-1/k)-1) * ( np.log(b+k*(iloss-160)) - np.log(b))
Last but not least, it's highly recommended to make yourself familiar with vectorized numpy operations to avoid loops:
def obj_function(x):
k, b = x
iloss = loss[loss > 160.0]
q = ((-1/k)-1)
target = np.log(1/b) + q* ( np.log(b+k*(iloss-160)) - np.log(b))
return -1.0*np.sum(target)
This gives me
In [9]: minimize(obj_function, x0=[0.7, 20])
Out[9]:
fun: 108.20609317144068
hess_inv: array([[ 1.03577269e-01, -2.40749570e+00],
[-2.40749570e+00, 1.45377397e+02]])
jac: array([ 9.53674316e-07, -9.53674316e-07])
message: 'Optimization terminated successfully.'
nfev: 39
nit: 9
njev: 13
status: 0
success: True
x: array([ 0.43624657, 32.53159127])