按照中给出的示例lmfit
,我正在尝试设置一个与我的问题类似的示例。我的问题最初是在我的数据中我可以拟合两个或三个模型,而我的模型是高度非线性的,但每个模型只有一个自由参数。
我的示例类似于lmfit
文档:
x = np.linspace(0, 15, 301)
data = (5. * np.sin(2 * x - 0.1) * np.exp(-x*x*0.025) +(-2.6 * np.sin(-0.6 * x + 1.5) * np.exp(-x*x*3.0)+np.random.normal(size=len(x), scale=0.2) ))
def fcn2min(params, x, data):
model=0
for i in range(2):
exec("amp_%d=%s"%(i+1,repr(params['amp_%d'%(i+1)].value)))
exec("shift_%d=%s"%(i+1,repr(params['shift_%d'%(i+1)].value)))
exec("omega_%d=%s"%(i+1,repr(params['omega_%d'%(i+1)].value)))
exec("decay_%d=%s"%(i+1,repr(params['decay_%d'%(i+1)].value)))
model += eval("amp_%d"%(i+1)) * np.sin(x * eval("omega_%d"%(i+1)) + eval("shift_%d"%(i+1))) * np.exp(-x*x*eval("decay_%d"%(i+1)))
return (model-data)/data
params=Parameters()
for i in range(2):
params.add('amp_%d'%(i+1), value= 10, vary=True, min=-3, max=3)
params.add('decay_%d'%(i+1), value= 0.1,vary=True,min=0,max=4.)
params.add('shift_%d'%(i+1), value= 0.0, vary=True,min=-np.pi, max=np.pi)
params.add('omega_%d'%(i+1), value= 3.0, vary=True,min=-2.5, max=2.5)
result = minimize(fcn2min, params, args=(x, data),method='nelder')
得到的结果:
final = data + result.residual
# write error report
report_fit(params)
[[Variables]]
amp_1: -1.74789852 (init= 3)
decay_1: 0.05493661 (init= 0.1)
shift_1: 0.07807319 (init= 0)
omega_1: -2.00291964 (init= 2.5)
amp_2: -1.30857699 (init= 3)
decay_2: 0.82303744 (init= 0.1)
shift_2: -0.04742474 (init= 0)
omega_2: 2.44085535 (init= 2.5)
[[Correlations]] (unreported correlations are < 0.100)
自由参数看起来完全不正确,但在最终结果图中,很明显它遵循数据的分布,但幅度并不完全正确
try:
import pylab
pylab.plot(x, data, 'k+')
pylab.plot(x, final, 'r')
pylab.show()
except:
pass
有什么建议可以修改代码以获得正确的结果吗?