我想用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’]
和的最佳值。所以你有了主意,接下来会发生什么
现在让我向您展示我是如何为此编写代码的。首先,我定义我的目标函数:
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文件的链接放在我解释的末尾,这样你可以自己检查一下
如下图所示,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。(您也可以在数据选项卡-->;分析-->;解决方案中看到优化选项。)
excel和python中用于优化的数据是相同的,可以使用here(非常简单,包含501行和1列)。
*如果你不能下载这些文件,请告诉我,然后我会更新它们
初始化将
[0.7, 20]
的值应用于每个参数。应使用标量来初始化value
,例如:另一个问题是
gekko
需要使用特殊函数来为解算器执行自动微分。对于目标函数,切换到求和的gekko
版本,如下所示:如果通过更改
loss
的值来计算x
,则目标函数具有logical expressions that need special treatment用于基于梯度的解算器的解。尝试对条件语句使用if3()
函数,或者slack variables(首选)。对目标函数求值一次,以构建符号表达式,然后将其编译为字节码,并使用其中一个解算器进行求解。符号表达式位于m.path
文件的gk0_model.apm
中对编辑的响应
感谢您发布完整代码的编辑。以下是一个潜在的解决方案:
解算器到达解的边界。边界可能需要加宽,以使任意限制不是解决方案
更新
这是最终的解决方案。代码中的目标函数有一个问题,因此应该进行修复。下面是正确的脚本:
输出:
结果与Microsoft Excel的优化结果非常接近
如果我能正确地看到,
My_DataFrame
已在全局范围中定义问题是
obj_funtion
尝试访问它(成功),然后修改它的值(失败) 这是因为默认情况下不能从局部范围修改全局变量修复:
在
obj_function
的开头添加一行:这会解决你的问题
附加说明:
如果您只想访问
My_DataFrame
,它将不会出现任何错误,并且您不需要添加global
关键字还有,我只是想感谢你为此付出的努力。这里有一个关于你想做什么的适当解释,相关的背景信息,一个优秀的图表(
Whiteboard
也非常好),甚至还有一个最小的工作示例。 所有这些问题都应该是这样的,这会让每个人的生活更轻松qw.Maximize()
只设置优化的目标,您仍然需要在模型上调用solve()
相关问题 更多 >
编程相关推荐