SciPy.optimize.least_squares() Objective 函数问题
SciPy.optimize.least_squares() Objective Function Questions
我试图通过优化三个未知参数 a、b 和 c0 来最小化高度非线性函数。我试图在 Python 3.
中复制赌场轮盘球的一些控制方程
这是研究论文的 link:
http://www.dewtronics.com/tutorials/roulette/documents/Roulette_Physik.pdf
我将在论文中引用等式 (35) 和 (40)。
基本上,我用秒表测量轮盘上旋转的轮盘球的圈数。对于每个连续的圈,由于非保守摩擦力的动量损失,圈时间将增加。然后,我使用等式 (40) 中的 Levenberg-Marquardt 最小二乘法进行这些时间测量并拟合等式 (35)。
我的问题有两个:
(1)我用的是scipy.optimize.least_squares()方法='lm',不知道objective函数怎么写!现在我已经完全按照论文中的方式编写了函数:
def fall_time(k,a,b,c0):
F = (1 / (a * b)) * (c0 - np.arcsinh(c0) * np.exp(a * k * 2 * np.pi))
return F
def parameter_estimation_function(x0,tk):
a = x0[0]
b = x0[1]
c0 = x0[2]
S = 0
for i,t in enumerate(tk):
k = i + 1
S += (t - fall_time(k,a,b,c0))**2
return [S,1,1]
sol = least_squares(parameter_estimation_function,[0.1,0.8,-0.1],args=([tk1]),method='lm',jac='2-point',max_nfev=2000)
print(sol)
现在,在文档示例中,我从未见过 objective 函数是按照我的方式编写的。在文档中,objective 函数始终是 returns 残差,而不是残差的平方。此外,在文档中他们从不使用总和!所以我想知道总和和平方是否在 least_squares()
的幕后自动处理?
(2) 也许我的第二个问题是因为我没有理解如何编写 objective 函数。但无论如何,我无法让算法收敛于最小值。我知道这是因为 levenberg 算法是 "greedy" 并在最接近的最小值附近停止,但我认为在不同的初始猜测下我至少能够收敛到大致相同的结果。对最初的猜测稍作改动,我得到了带有不同符号的参数结果。此外,我还没有找到允许算法收敛的初始猜测组合!它总是在找到解决方案之前超时。我什至将函数评估的数量增加到 10,000,看看它是否会。无济于事!
也许有人可以在这里阐明我的错误!我对 python 和 scipy 库还比较陌生!
这是我从此处的视频中自己测量的 tk
的一些样本数据:https://www.youtube.com/watch?v=0Zj_9ypBnzg
tk = [0.52,1.28,2.04,3.17,4.53,6.22]
tk1 = [0.51,1.4,2.09,3,4.42,6.17]
tk2 = [0.63,1.35,2.19,3.02,4.57,6.29]
tk3 = [0.63,1.39,2.23,3.28,4.70,6.32]
tk4 = [0.57,1.4,2.1,3.06,4.53,6.17]
谢谢
1) 是的,正如您所怀疑的,残差的总和和平方是自动处理的。
2) 很难说,因为我对这个问题不是很熟悉(例如,存在多少个局部最小值,什么构成 'reasonable' 结果等)。稍后我可能会调查更多。
但为了好玩,我摆弄了一些值,看看会发生什么。例如,您可以将 1/b
常量替换为独立变量 b_inv
,这似乎使结果稳定了很多。这是我用来检查结果的代码。 (请注意,为简洁起见,我重写了 objective 函数。它只是利用了 numpy 数组的逐元素操作,而没有改变整体结果。)
import numpy as np
from scipy.optimize import least_squares
def fall_time(k,a,b_inv,c0):
return (b_inv / a) * (c0 - np.arcsinh(c0) * np.exp(a * k * 2 * np.pi))
def parameter_estimation_function(x,tk):
return np.asarray(tk) - fall_time(k=np.arange(1,len(tk)+1), a=x[0],b_inv=x[1],c0=x[2])
tk_samples = [
[0.52,1.28,2.04,3.17,4.53,6.22],
[0.51,1.4,2.09,3,4.42,6.17],
[0.63,1.35,2.19,3.02,4.57,6.29],
[0.63,1.39,2.23,3.28,4.70,6.32],
[0.57,1.4,2.1,3.06,4.53,6.17]
]
for i in range(len(tk_samples)):
sol = least_squares(parameter_estimation_function,[0.1,1.25,-0.1],
args=(tk_samples[i],),method='lm',jac='2-point',max_nfev=2000)
print(sol.x)
控制台输出:
[ 0.03621789 0.64201913 -0.12072879]
[ 3.59319972e-02 1.17129458e+01 -6.53358716e-03]
[ 3.55516005e-02 1.48491493e+01 -5.31098257e-03]
[ 3.18068316e-02 1.11828091e+01 -7.75329834e-03]
[ 3.43920725e-02 1.25160378e+01 -6.36307506e-03]
我试图通过优化三个未知参数 a、b 和 c0 来最小化高度非线性函数。我试图在 Python 3.
中复制赌场轮盘球的一些控制方程这是研究论文的 link: http://www.dewtronics.com/tutorials/roulette/documents/Roulette_Physik.pdf
我将在论文中引用等式 (35) 和 (40)。 基本上,我用秒表测量轮盘上旋转的轮盘球的圈数。对于每个连续的圈,由于非保守摩擦力的动量损失,圈时间将增加。然后,我使用等式 (40) 中的 Levenberg-Marquardt 最小二乘法进行这些时间测量并拟合等式 (35)。
我的问题有两个: (1)我用的是scipy.optimize.least_squares()方法='lm',不知道objective函数怎么写!现在我已经完全按照论文中的方式编写了函数:
def fall_time(k,a,b,c0):
F = (1 / (a * b)) * (c0 - np.arcsinh(c0) * np.exp(a * k * 2 * np.pi))
return F
def parameter_estimation_function(x0,tk):
a = x0[0]
b = x0[1]
c0 = x0[2]
S = 0
for i,t in enumerate(tk):
k = i + 1
S += (t - fall_time(k,a,b,c0))**2
return [S,1,1]
sol = least_squares(parameter_estimation_function,[0.1,0.8,-0.1],args=([tk1]),method='lm',jac='2-point',max_nfev=2000)
print(sol)
现在,在文档示例中,我从未见过 objective 函数是按照我的方式编写的。在文档中,objective 函数始终是 returns 残差,而不是残差的平方。此外,在文档中他们从不使用总和!所以我想知道总和和平方是否在 least_squares()
的幕后自动处理?
(2) 也许我的第二个问题是因为我没有理解如何编写 objective 函数。但无论如何,我无法让算法收敛于最小值。我知道这是因为 levenberg 算法是 "greedy" 并在最接近的最小值附近停止,但我认为在不同的初始猜测下我至少能够收敛到大致相同的结果。对最初的猜测稍作改动,我得到了带有不同符号的参数结果。此外,我还没有找到允许算法收敛的初始猜测组合!它总是在找到解决方案之前超时。我什至将函数评估的数量增加到 10,000,看看它是否会。无济于事!
也许有人可以在这里阐明我的错误!我对 python 和 scipy 库还比较陌生!
这是我从此处的视频中自己测量的 tk
的一些样本数据:https://www.youtube.com/watch?v=0Zj_9ypBnzg
tk = [0.52,1.28,2.04,3.17,4.53,6.22]
tk1 = [0.51,1.4,2.09,3,4.42,6.17]
tk2 = [0.63,1.35,2.19,3.02,4.57,6.29]
tk3 = [0.63,1.39,2.23,3.28,4.70,6.32]
tk4 = [0.57,1.4,2.1,3.06,4.53,6.17]
谢谢
1) 是的,正如您所怀疑的,残差的总和和平方是自动处理的。
2) 很难说,因为我对这个问题不是很熟悉(例如,存在多少个局部最小值,什么构成 'reasonable' 结果等)。稍后我可能会调查更多。
但为了好玩,我摆弄了一些值,看看会发生什么。例如,您可以将 1/b
常量替换为独立变量 b_inv
,这似乎使结果稳定了很多。这是我用来检查结果的代码。 (请注意,为简洁起见,我重写了 objective 函数。它只是利用了 numpy 数组的逐元素操作,而没有改变整体结果。)
import numpy as np
from scipy.optimize import least_squares
def fall_time(k,a,b_inv,c0):
return (b_inv / a) * (c0 - np.arcsinh(c0) * np.exp(a * k * 2 * np.pi))
def parameter_estimation_function(x,tk):
return np.asarray(tk) - fall_time(k=np.arange(1,len(tk)+1), a=x[0],b_inv=x[1],c0=x[2])
tk_samples = [
[0.52,1.28,2.04,3.17,4.53,6.22],
[0.51,1.4,2.09,3,4.42,6.17],
[0.63,1.35,2.19,3.02,4.57,6.29],
[0.63,1.39,2.23,3.28,4.70,6.32],
[0.57,1.4,2.1,3.06,4.53,6.17]
]
for i in range(len(tk_samples)):
sol = least_squares(parameter_estimation_function,[0.1,1.25,-0.1],
args=(tk_samples[i],),method='lm',jac='2-point',max_nfev=2000)
print(sol.x)
控制台输出:
[ 0.03621789 0.64201913 -0.12072879]
[ 3.59319972e-02 1.17129458e+01 -6.53358716e-03]
[ 3.55516005e-02 1.48491493e+01 -5.31098257e-03]
[ 3.18068316e-02 1.11828091e+01 -7.75329834e-03]
[ 3.43920725e-02 1.25160378e+01 -6.36307506e-03]