python curve_fit 不适用于刚性模型
python curve_fit doesn't work with stiff model
我正在尝试找到尽可能适合绿色曲线上的蓝色模型的x0
参数(x0
控制垛口的宽度;见下文)。
这是我的尝试:
from pylab import *
from scipy.optimize import curve_fit
x=linspace(0,2*pi,1000)
def crenel(x):return sign(sin(x))
def inverter(x,x0): return (crenel(x-x0)+crenel(x+x0))/2
p,e = curve_fit(inverter,x,sin(x),1)
plot(x,inverter(x,*p),x,sin(x))
ylim(-1.5,1.5)
手算,最优值为 x0 = arcsin(1/2) # 0.523598
,但 curve_fit 未估计任何值 ("OptimizeWarning: Covariance of the parameters could not be estimated")。我怀疑模型的刚度。 docs 告知:
The algorithm uses the Levenberg-Marquardt algorithm through leastsq. Additional keyword arguments are passed directly to that algorithm.
所以我的问题是:在这种情况下,是否有关键字参数可以帮助 curve_fit
估计参数?或另一种方法?
感谢您的任何建议。
问题是 curve_fit
试图最小化的 objective 函数不连续。 x0
控制 inverter
函数中不连续点的位置。当不连续点穿过 x
中的一个网格点时,objective 函数中会发生跳跃。在这些点之间,objective 函数是常量。 curve_fit
(其实是leastsq
,curve_fit
使用的函数)并不是为了处理这样的函数而设计的
以下函数 sse
是(实际上)curve_fit
试图最小化的函数,x
与示例中定义的 x
相同,并且y = sin(x)
:
def sse(x0, x, y):
f = inverter(x, x0)
diff = y - f
s = (diff**2).sum()
return s
如果您使用
等代码将此函数绘制在精细网格上
xx = np.linspace(0, 1, 10000)
yy = [sse(x0, x, y) for x0 in xx]
plot(xx, yy)
然后放大,你会看到
要使用 scipy
找到最佳值,您可以使用 fmin
和平滑的 objective 函数。例如,这里是连续 objective 函数,仅使用区间 [0, pi/2](quad
是 scipy.integrate.quad
):
def func(x0):
s0, e0 = quad(lambda x: np.sin(x)**2, 0, x0)
s1, e0 = quad(lambda x: (1 - np.sin(x))**2, x0, 0.5*np.pi)
return s0 + s1
scipy.optimize.fmin
可用于查找该函数的最小值,如 ipython 会话中的此片段所示:
In [202]: fmin(func, 0.3, xtol=1e-8)
Optimization terminated successfully.
Current function value: 0.100545
Iterations: 28
Function evaluations: 56
Out[202]: array([ 0.52359878])
In [203]: np.arcsin(0.5)
Out[203]: 0.52359877559829882
我正在尝试找到尽可能适合绿色曲线上的蓝色模型的x0
参数(x0
控制垛口的宽度;见下文)。
这是我的尝试:
from pylab import *
from scipy.optimize import curve_fit
x=linspace(0,2*pi,1000)
def crenel(x):return sign(sin(x))
def inverter(x,x0): return (crenel(x-x0)+crenel(x+x0))/2
p,e = curve_fit(inverter,x,sin(x),1)
plot(x,inverter(x,*p),x,sin(x))
ylim(-1.5,1.5)
手算,最优值为 x0 = arcsin(1/2) # 0.523598
,但 curve_fit 未估计任何值 ("OptimizeWarning: Covariance of the parameters could not be estimated")。我怀疑模型的刚度。 docs 告知:
The algorithm uses the Levenberg-Marquardt algorithm through leastsq. Additional keyword arguments are passed directly to that algorithm.
所以我的问题是:在这种情况下,是否有关键字参数可以帮助 curve_fit
估计参数?或另一种方法?
感谢您的任何建议。
问题是 curve_fit
试图最小化的 objective 函数不连续。 x0
控制 inverter
函数中不连续点的位置。当不连续点穿过 x
中的一个网格点时,objective 函数中会发生跳跃。在这些点之间,objective 函数是常量。 curve_fit
(其实是leastsq
,curve_fit
使用的函数)并不是为了处理这样的函数而设计的
以下函数 sse
是(实际上)curve_fit
试图最小化的函数,x
与示例中定义的 x
相同,并且y = sin(x)
:
def sse(x0, x, y):
f = inverter(x, x0)
diff = y - f
s = (diff**2).sum()
return s
如果您使用
等代码将此函数绘制在精细网格上xx = np.linspace(0, 1, 10000)
yy = [sse(x0, x, y) for x0 in xx]
plot(xx, yy)
然后放大,你会看到
要使用 scipy
找到最佳值,您可以使用 fmin
和平滑的 objective 函数。例如,这里是连续 objective 函数,仅使用区间 [0, pi/2](quad
是 scipy.integrate.quad
):
def func(x0):
s0, e0 = quad(lambda x: np.sin(x)**2, 0, x0)
s1, e0 = quad(lambda x: (1 - np.sin(x))**2, x0, 0.5*np.pi)
return s0 + s1
scipy.optimize.fmin
可用于查找该函数的最小值,如 ipython 会话中的此片段所示:
In [202]: fmin(func, 0.3, xtol=1e-8)
Optimization terminated successfully.
Current function value: 0.100545
Iterations: 28
Function evaluations: 56
Out[202]: array([ 0.52359878])
In [203]: np.arcsin(0.5)
Out[203]: 0.52359877559829882