python 中 scipy.optimize 的 MLE 应用程序

MLE application with scipy.optimize in python

我想运行在python中进行简单的最大似然估计。我想在 python 中使用 Scipy.optimize.minimize 来尝试一下。首先,我将解释我的模型,以便您了解将要发生的事情。

模型说明

通过 MLE 我想估计最大化我的 objective 函数的 2 个变量的最佳值!
这些变量是什么以及 objective 函数是什么?在下面的图片中你可以看到我的 objective 功能!我的变量是 'Betaa''Kasi'.

现在你知道我的 objective 函数是什么,我的变量是什么(objective 函数中的 item 除外。我会解释它!请继续阅读)。我想在数据集上实现 MLE(click 从我的 google 驱动器下载它),你可以在下面的图片中看到:

这是非常简单的数据集!包含 501 行和 1 列。现在您已经看到了我的数据,是时候解释 MLE 如何处理它了!

MLE

对于 'Loss' 列中的每个值,如果该值 >160,我应该通过 objective 函数计算目标值(我在第一部分中谈到过),如果该值是<160,目标为 0!看下面的图片然后你会发现我在说什么(你也会发现我说的 item 是什么)! MLE 将通过为 'Betaa' 和 'Kasi'!

找到最佳值来最大化 Target 中的值总和

代码

所以这是为此目的编写的脚本:

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)

此代码的结果:

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

正如您在结果中看到的那样,success: False 并且未找到最佳解决方案!
现在我不知道如何解决这个问题以及如何解决这个 MLE 问题。
任何帮助将不胜感激。

正如评论中已经提到的那样,在评估功率时由于溢出而出现多个运行时警告。仔细观察,您会发现您的 objective 函数是错误的。根据你的图片,应该是

(1 + ( k*(iloss-160)/b) )**( (-1/k) - 1 ) 

而不是

(1 + ( k*(iloss-160)/b )**( (-1/k) - 1 )

此外,您可以避免计算幂并使用对数规则重写您的objective:

np.log(1/b) + ((-1/k)-1) * ( np.log(b+k*(iloss-160)) - np.log(b))

最后但同样重要的是,强烈建议让自己熟悉矢量化 numpy 操作以避免循环:

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)

这给了我

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