python 中带有 gekko 的 MLE 应用程序
MLE application with gekko in python
我想在 python 中使用 gekko 包实现 MLE（最大似然估计）。假设我们有一个包含两列的 DataFrame：['Loss'、'Target'] 并且它的长度等于 500。
首先我们必须导入我们需要的包：
from gekko import GEKKO
import numpy as np
import pandas as pd
然后我们简单地创建 DataFrame 这样的：
My_DataFrame = pd.DataFrame({"Loss":np.linspace(-555.795 , 477.841 , 500) , "Target":0.0})
My_DataFrame = My_DataFrame.sort_values(by=["Loss"] , ascending=False).reset_index(drop=True)
My_DataFrame
它看起来像这样：
['Target']列的一些组成部分应该用我在图片下方写下的公式计算（其余部分保持为零。我继续解释更多，请继续阅读）这样你就可以完美地看到它。配方的两个主要元素是“Kasi”和“Betaa”。我想为他们找到最大化 My_DataFrame[‘Target’] 之和的最佳价值。所以你明白了会发生什么！
现在让我向您展示我是如何为此目的编写代码的。首先我定义我的 objective 函数：
def obj_function(Array):
"""
[Purpose]:
+ it will calculate each component of My_DataFrame["Target"] column! then i can maximize sum(My_DataFrame["Target"]) and find best 'Kasi' and 'Betaa' for it!
[Parameters]:
+ This function gets Array that contains 'Kasi' and 'Betaa'.
Array[0] represents 'Kasi' and Array[1] represents 'Betaa'
[returns]:
+ returns a pandas.series.
actually it returns new components of My_DataFrame["Target"]
"""
# in following code if you don't know what is `qw`, just look at the next code cell right after this cell (I mean next section).
# in following code np.where(My_DataFrame["Loss"] == item)[0][0] is telling me the row's index of item.
for item in My_DataFrame[My_DataFrame["Loss"]>160]['Loss']:
My_DataFrame.iloc[np.where(My_DataFrame["Loss"] == item)[0][0] , 1] = qw.log10((1/Array[1])*( 1 + (Array[0]*(item-160)/Array[1])**( (-1/Array[0]) - 1 )))
return My_DataFrame["Target"]
如果您对 obj_function 函数中的 for loop 发生了什么感到困惑，请查看下图，其中包含一个简短示例！如果没有，请跳过这部分：
那么我们只需要进行优化即可。为此，我使用 gekko 包。 请注意我想找到“Kasi”和“Betaa”的最佳值，所以我有两个主要变量而且我没有任何限制！
那么让我们开始吧：
# i have 2 variables : 'Kasi' and 'Betaa', so I put nd=2
nd = 2
qw = GEKKO()
# now i want to specify my variables ('Kasi' and 'Betaa') with initial values --> Kasi = 0.7 and Betaa = 20.0
x = qw.Array(qw.Var , nd , value = [0.7 , 20])
# So i guess now x[0] represents 'Kasi' and x[1] represents 'Betaa'
qw.Maximize(np.sum(obj_function(x)))
然后当我想用qw.solve()解决优化时：
qw.solve()
但是我得到了这个错误：
Exception: This steady-state IMODE only allows scalar values.
我该如何解决这个问题？ （为了方便起见，下一节收集了完整的脚本）
from gekko import GEKKO
import numpy as np
import pandas as pd
My_DataFrame = pd.DataFrame({"Loss":np.linspace(-555.795 , 477.841 , 500) , "Target":0.0})
My_DataFrame = My_DataFrame.sort_values(by=["Loss"] , ascending=False).reset_index(drop=True)
def obj_function(Array):
"""
[Purpose]:
+ it will calculate each component of My_DataFrame["Target"] column! then i can maximize sum(My_DataFrame["Target"]) and find best 'Kasi' and 'Betaa' for it!
[Parameters]:
+ This function gets Array that contains 'Kasi' and 'Betaa'.
Array[0] represents 'Kasi' and Array[1] represents 'Betaa'
[returns]:
+ returns a pandas.series.
actually it returns new components of My_DataFrame["Target"]
"""
# in following code if you don't know what is `qw`, just look at the next code cell right after this cell (I mean next section).
# in following code np.where(My_DataFrame["Loss"] == item)[0][0] is telling me the row's index of item.
for item in My_DataFrame[My_DataFrame["Loss"]>160]['Loss']:
My_DataFrame.iloc[np.where(My_DataFrame["Loss"] == item)[0][0] , 1] = qw.log10((1/Array[1])*( 1 + (Array[0]*(item-160)/Array[1])**( (-1/Array[0]) - 1 )))
return My_DataFrame["Target"]
# i have 2 variables : 'Kasi' and 'Betaa', so I put nd=2
nd = 2
qw = GEKKO()
# now i want to specify my variables ('Kasi' and 'Betaa') with initial values --> Kasi = 0.7 and Betaa = 20.0
x = qw.Array(qw.Var , nd)
for i,xi in enumerate([0.7, 20]):
x[i].value = xi
# So i guess now x[0] represents 'Kasi' and x[1] represents 'Betaa'
qw.Maximize(qw.sum(obj_function(x)))
提议的潜在脚本在这里：
from gekko import GEKKO
import numpy as np
import pandas as pd
My_DataFrame = pd.read_excel("[<FILE_PATH_IN_YOUR_MACHINE>]\Losses.xlsx")
# i'll put link of "Losses.xlsx" file in the end of my explaination
# so you can download it from my google drive.
loss = My_DataFrame["Loss"]
def obj_function(x):
k,b = x
target = []
for iloss in loss:
if iloss>160:
t = qw.log((1/b)*(1+(k*(iloss-160)/b)**((-1/k)-1)))
target.append(t)
return target
qw = GEKKO(remote=False)
nd = 2
x = qw.Array(qw.Var,nd)
# initial values --> Kasi = 0.7 and Betaa = 20.0
for i,xi in enumerate([0.7, 20]):
x[i].value = xi
# bounds
k,b = x
k.lower=0.1; k.upper=0.8
b.lower=10; b.upper=500
qw.Maximize(qw.sum(obj_function(x)))
qw.options.SOLVER = 1
qw.solve()
print('k = ',k.value[0])
print('b = ',b.value[0])
python 输出：
objective function = -1155.4861315885942
b = 500.0
k = 0.1
请注意 python 输出中的 b 代表“Betaa”，k 代表“Kasi”。
输出看起来有点奇怪，所以 我决定测试一下！ 为此我使用了 Microsoft Excel Solver!
（我把 excel 文件的 link 放在我解释的最后，所以你可以自己检查一下
你想要的。）正如你在下面的图片中看到的，excel 的优化已经完成并且最优解
已成功找到（优化见下图结果）。
excel 输出：
objective function = -108.21
Betaa = 32.53161
Kasi = 0.436246
如您所见，python output 和 excel output 之间存在巨大差异，似乎 excel 表现得相当不错！ 所以我想问题仍然存在，建议 python 脚本表现不佳...
Implementation_in_Excel.xls Microsoft excel 应用程序的优化文件可用 here。（您还可以在数据选项卡 --> 分析 --> Slover 中查看优化选项。）
excel 和 python 中用于优化的数据相同并且可用 here （非常简单，包含 501 行和 1 列）。
*如果你不能下载文件，让我知道，我会更新它们。
如果我没看错的话，My_DataFrame已经在全局范围内定义了。
问题是 obj_funtion 尝试访问它（成功）然后修改它的值（失败）
这是因为默认情况下您无法从局部范围修改全局变量。
修复：
在obj_function开头添加一行：
def obj_function(Array):
# comments
global My_DataFrame
for item .... # remains same
这应该可以解决您的问题。
补充说明：
如果您只想访问 My_DataFrame，它可以正常工作，并且您不需要添加 global 关键字
此外，我只是想感谢您为此付出的努力。有对您想做什么的正确解释、相关的背景信息、出色的图表（Whiteboard 也非常棒），甚至还有一个最小的工作示例。
这应该是所有SO问题的方式，它会让每个人的生活更轻松
qw.Maximize()只设置优化的objective，你还需要在你的模型上调用solve()。
初始化正在将 [0.7, 20] 的值应用于每个参数。应该使用标量来初始化 value 而不是：
x = qw.Array(qw.Var , nd)
for i,xi in enumerate([0.7, 20]):
x[i].value = xi
另一个问题是gekko需要使用特殊函数对求解器进行自动微分。对于 objective 函数，切换到求和的 gekko 版本：
qw.Maximize(qw.sum(obj_function(x)))
如果 loss 是通过更改 x 的值计算的，则 objective 函数具有 logical expressions that need special treatment for solution with gradient-based solvers. Try using the if3() function for a conditional statement or else slack variables（首选）。 objective 函数被评估一次以构建一个符号表达式，然后将其编译为字节代码并使用其中一个求解器求解。符号表达式位于 gk0_model.apm 文件的 m.path 中。
回复编辑
感谢您发布包含完整代码的编辑。这是一个潜在的解决方案：
from gekko import GEKKO
import numpy as np
import pandas as pd
loss = np.linspace(-555.795 , 477.841 , 500)
def obj_function(x):
k,b = x
target = []
for iloss in loss:
if iloss>160:
t = qw.log((1/b)*(1+(k*(iloss-160)/b)**((-1/k)-1)))
target.append(t)
return target
qw = GEKKO(remote=False)
nd = 2
x = qw.Array(qw.Var,nd)
# initial values --> Kasi = 0.7 and Betaa = 20.0
for i,xi in enumerate([0.7, 20]):
x[i].value = xi
# bounds
k,b = x
k.lower=0.6; k.upper=0.8
b.lower=10; b.upper=30
qw.Maximize(qw.sum(obj_function(x)))
qw.options.SOLVER = 1
qw.solve()
print('k = ',k.value[0])
print('b = ',b.value[0])
求解器到达解的边界。可能需要扩大范围，这样任意限制就不是解决方案。
更新
这是最终的解决方案。代码中的 objective 函数有问题所以应该修复这是正确的脚本：
from gekko import GEKKO
import numpy as np
import pandas as pd
My_DataFrame = pd.read_excel("<FILE_PATH_IN_YOUR_MACHINE>\Losses.xlsx")
loss = My_DataFrame["Loss"]
def obj_function(x):
k,b = x
q = ((-1/k)-1)
target = []
for iloss in loss:
if iloss>160:
t = qw.log(1/b) + q* ( qw.log(b+k*(iloss-160)) - qw.log(b))
target.append(t)
return target
qw = GEKKO(remote=False)
nd = 2
x = qw.Array(qw.Var,nd)
# initial values --> Kasi = 0.7 and Betaa = 20.0
for i,xi in enumerate([0.7, 20]):
x[i].value = xi
qw.Maximize(qw.sum(obj_function(x)))
qw.solve()
print('Kasi = ',x[0].value)
print('Betaa = ',x[1].value)
输出：
The final value of the objective function is 108.20609317143486
---------------------------------------------------
Solver : IPOPT (v3.12)
Solution time : 0.031200000000000006 sec
Objective : 108.20609317143486
Successful solution
---------------------------------------------------
Kasi = [0.436245842]
Betaa = [32.531632983]
结果接近微软的优化结果Excel。
我想在 python 中使用 gekko 包实现 MLE（最大似然估计）。假设我们有一个包含两列的 DataFrame：['Loss'、'Target'] 并且它的长度等于 500。
首先我们必须导入我们需要的包：
from gekko import GEKKO
import numpy as np
import pandas as pd
然后我们简单地创建 DataFrame 这样的：
My_DataFrame = pd.DataFrame({"Loss":np.linspace(-555.795 , 477.841 , 500) , "Target":0.0})
My_DataFrame = My_DataFrame.sort_values(by=["Loss"] , ascending=False).reset_index(drop=True)
My_DataFrame
它看起来像这样：
['Target']列的一些组成部分应该用我在图片下方写下的公式计算（其余部分保持为零。我继续解释更多，请继续阅读）这样你就可以完美地看到它。配方的两个主要元素是“Kasi”和“Betaa”。我想为他们找到最大化 My_DataFrame[‘Target’] 之和的最佳价值。所以你明白了会发生什么！
现在让我向您展示我是如何为此目的编写代码的。首先我定义我的 objective 函数：
def obj_function(Array):
"""
[Purpose]:
+ it will calculate each component of My_DataFrame["Target"] column! then i can maximize sum(My_DataFrame["Target"]) and find best 'Kasi' and 'Betaa' for it!
[Parameters]:
+ This function gets Array that contains 'Kasi' and 'Betaa'.
Array[0] represents 'Kasi' and Array[1] represents 'Betaa'
[returns]:
+ returns a pandas.series.
actually it returns new components of My_DataFrame["Target"]
"""
# in following code if you don't know what is `qw`, just look at the next code cell right after this cell (I mean next section).
# in following code np.where(My_DataFrame["Loss"] == item)[0][0] is telling me the row's index of item.
for item in My_DataFrame[My_DataFrame["Loss"]>160]['Loss']:
My_DataFrame.iloc[np.where(My_DataFrame["Loss"] == item)[0][0] , 1] = qw.log10((1/Array[1])*( 1 + (Array[0]*(item-160)/Array[1])**( (-1/Array[0]) - 1 )))
return My_DataFrame["Target"]
如果您对 obj_function 函数中的 for loop 发生了什么感到困惑，请查看下图，其中包含一个简短示例！如果没有，请跳过这部分：
那么我们只需要进行优化即可。为此，我使用 gekko 包。 请注意我想找到“Kasi”和“Betaa”的最佳值，所以我有两个主要变量而且我没有任何限制！
那么让我们开始吧：
# i have 2 variables : 'Kasi' and 'Betaa', so I put nd=2
nd = 2
qw = GEKKO()
# now i want to specify my variables ('Kasi' and 'Betaa') with initial values --> Kasi = 0.7 and Betaa = 20.0
x = qw.Array(qw.Var , nd , value = [0.7 , 20])
# So i guess now x[0] represents 'Kasi' and x[1] represents 'Betaa'
qw.Maximize(np.sum(obj_function(x)))
然后当我想用qw.solve()解决优化时：
qw.solve()
但是我得到了这个错误：
Exception: This steady-state IMODE only allows scalar values.
我该如何解决这个问题？ （为了方便起见，下一节收集了完整的脚本）
from gekko import GEKKO
import numpy as np
import pandas as pd
My_DataFrame = pd.DataFrame({"Loss":np.linspace(-555.795 , 477.841 , 500) , "Target":0.0})
My_DataFrame = My_DataFrame.sort_values(by=["Loss"] , ascending=False).reset_index(drop=True)
def obj_function(Array):
"""
[Purpose]:
+ it will calculate each component of My_DataFrame["Target"] column! then i can maximize sum(My_DataFrame["Target"]) and find best 'Kasi' and 'Betaa' for it!
[Parameters]:
+ This function gets Array that contains 'Kasi' and 'Betaa'.
Array[0] represents 'Kasi' and Array[1] represents 'Betaa'
[returns]:
+ returns a pandas.series.
actually it returns new components of My_DataFrame["Target"]
"""
# in following code if you don't know what is `qw`, just look at the next code cell right after this cell (I mean next section).
# in following code np.where(My_DataFrame["Loss"] == item)[0][0] is telling me the row's index of item.
for item in My_DataFrame[My_DataFrame["Loss"]>160]['Loss']:
My_DataFrame.iloc[np.where(My_DataFrame["Loss"] == item)[0][0] , 1] = qw.log10((1/Array[1])*( 1 + (Array[0]*(item-160)/Array[1])**( (-1/Array[0]) - 1 )))
return My_DataFrame["Target"]
# i have 2 variables : 'Kasi' and 'Betaa', so I put nd=2
nd = 2
qw = GEKKO()
# now i want to specify my variables ('Kasi' and 'Betaa') with initial values --> Kasi = 0.7 and Betaa = 20.0
x = qw.Array(qw.Var , nd)
for i,xi in enumerate([0.7, 20]):
x[i].value = xi
# So i guess now x[0] represents 'Kasi' and x[1] represents 'Betaa'
qw.Maximize(qw.sum(obj_function(x)))
提议的潜在脚本在这里：
from gekko import GEKKO
import numpy as np
import pandas as pd
My_DataFrame = pd.read_excel("[<FILE_PATH_IN_YOUR_MACHINE>]\Losses.xlsx")
# i'll put link of "Losses.xlsx" file in the end of my explaination
# so you can download it from my google drive.
loss = My_DataFrame["Loss"]
def obj_function(x):
k,b = x
target = []
for iloss in loss:
if iloss>160:
t = qw.log((1/b)*(1+(k*(iloss-160)/b)**((-1/k)-1)))
target.append(t)
return target
qw = GEKKO(remote=False)
nd = 2
x = qw.Array(qw.Var,nd)
# initial values --> Kasi = 0.7 and Betaa = 20.0
for i,xi in enumerate([0.7, 20]):
x[i].value = xi
# bounds
k,b = x
k.lower=0.1; k.upper=0.8
b.lower=10; b.upper=500
qw.Maximize(qw.sum(obj_function(x)))
qw.options.SOLVER = 1
qw.solve()
print('k = ',k.value[0])
print('b = ',b.value[0])
python 输出：
objective function = -1155.4861315885942
b = 500.0
k = 0.1
请注意 python 输出中的 b 代表“Betaa”，k 代表“Kasi”。
输出看起来有点奇怪，所以 我决定测试一下！ 为此我使用了 Microsoft Excel Solver!
（我把 excel 文件的 link 放在我解释的最后，所以你可以自己检查一下
你想要的。）正如你在下面的图片中看到的，excel 的优化已经完成并且最优解
已成功找到（优化见下图结果）。
excel 输出：
objective function = -108.21
Betaa = 32.53161
Kasi = 0.436246
如您所见，python output 和 excel output 之间存在巨大差异，似乎 excel 表现得相当不错！ 所以我想问题仍然存在，建议 python 脚本表现不佳...
Implementation_in_Excel.xls Microsoft excel 应用程序的优化文件可用 here。（您还可以在数据选项卡 --> 分析 --> Slover 中查看优化选项。）
excel 和 python 中用于优化的数据相同并且可用 here （非常简单，包含 501 行和 1 列）。
*如果你不能下载文件，让我知道，我会更新它们。
如果我没看错的话，My_DataFrame已经在全局范围内定义了。
问题是 obj_funtion 尝试访问它（成功）然后修改它的值（失败）
这是因为默认情况下您无法从局部范围修改全局变量。
修复：
在obj_function开头添加一行：
def obj_function(Array):
# comments
global My_DataFrame
for item .... # remains same
这应该可以解决您的问题。
补充说明：
如果您只想访问 My_DataFrame，它可以正常工作，并且您不需要添加 global 关键字
此外，我只是想感谢您为此付出的努力。有对您想做什么的正确解释、相关的背景信息、出色的图表（Whiteboard 也非常棒），甚至还有一个最小的工作示例。
这应该是所有SO问题的方式，它会让每个人的生活更轻松
qw.Maximize()只设置优化的objective，你还需要在你的模型上调用solve()。
初始化正在将 [0.7, 20] 的值应用于每个参数。应该使用标量来初始化 value 而不是：
x = qw.Array(qw.Var , nd)
for i,xi in enumerate([0.7, 20]):
x[i].value = xi
另一个问题是gekko需要使用特殊函数对求解器进行自动微分。对于 objective 函数，切换到求和的 gekko 版本：
qw.Maximize(qw.sum(obj_function(x)))
如果 loss 是通过更改 x 的值计算的，则 objective 函数具有 logical expressions that need special treatment for solution with gradient-based solvers. Try using the if3() function for a conditional statement or else slack variables（首选）。 objective 函数被评估一次以构建一个符号表达式，然后将其编译为字节代码并使用其中一个求解器求解。符号表达式位于 gk0_model.apm 文件的 m.path 中。
回复编辑
感谢您发布包含完整代码的编辑。这是一个潜在的解决方案：
from gekko import GEKKO
import numpy as np
import pandas as pd
loss = np.linspace(-555.795 , 477.841 , 500)
def obj_function(x):
k,b = x
target = []
for iloss in loss:
if iloss>160:
t = qw.log((1/b)*(1+(k*(iloss-160)/b)**((-1/k)-1)))
target.append(t)
return target
qw = GEKKO(remote=False)
nd = 2
x = qw.Array(qw.Var,nd)
# initial values --> Kasi = 0.7 and Betaa = 20.0
for i,xi in enumerate([0.7, 20]):
x[i].value = xi
# bounds
k,b = x
k.lower=0.6; k.upper=0.8
b.lower=10; b.upper=30
qw.Maximize(qw.sum(obj_function(x)))
qw.options.SOLVER = 1
qw.solve()
print('k = ',k.value[0])
print('b = ',b.value[0])
求解器到达解的边界。可能需要扩大范围，这样任意限制就不是解决方案。
更新
这是最终的解决方案。代码中的 objective 函数有问题所以应该修复这是正确的脚本：
from gekko import GEKKO
import numpy as np
import pandas as pd
My_DataFrame = pd.read_excel("<FILE_PATH_IN_YOUR_MACHINE>\Losses.xlsx")
loss = My_DataFrame["Loss"]
def obj_function(x):
k,b = x
q = ((-1/k)-1)
target = []
for iloss in loss:
if iloss>160:
t = qw.log(1/b) + q* ( qw.log(b+k*(iloss-160)) - qw.log(b))
target.append(t)
return target
qw = GEKKO(remote=False)
nd = 2
x = qw.Array(qw.Var,nd)
# initial values --> Kasi = 0.7 and Betaa = 20.0
for i,xi in enumerate([0.7, 20]):
x[i].value = xi
qw.Maximize(qw.sum(obj_function(x)))
qw.solve()
print('Kasi = ',x[0].value)
print('Betaa = ',x[1].value)
输出：
The final value of the objective function is 108.20609317143486
---------------------------------------------------
Solver : IPOPT (v3.12)
Solution time : 0.031200000000000006 sec
Objective : 108.20609317143486
Successful solution
---------------------------------------------------
Kasi = [0.436245842]
Betaa = [32.531632983]
结果接近微软的优化结果Excel。