pythonnumpyscipymaximizemle

MLE application with scipy.optimize in python


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.

Model Explanation

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'.
enter image description here
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:
enter image description here
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!

MLE

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)! enter image description here MLE is gonna maximize summation of values in Target by finding best values for 'Betaa' and 'Kasi'!

Code

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.


Solution

  • 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])