Python中的模拟退火拟合

Question

我正在尝试熟悉非线性拟合过程：双重退火。为此，我生成了一些合成数据，并尝试为它们拟合一个基本的 Furth 公式，请参见下面的代码：

import numpy as np
from numpy import savetxt
from numpy import genfromtxt
import random
import matplotlib.pyplot as plt 
from scipy.optimize import curve_fit
from lmfit import Minimizer, Parameters, report_fit
from scipy.optimize import dual_annealing

def Furth(A,B, a, x):
    model = A*x + B * (np.exp(-a*x)-1) 
    return model 

生成合成数据：

x = np.arange(200)*0.5
x = x[1:]
A =1.88
B = 2.35
a = 5602
y = Furth(A,B,a,x) + np.random.randn(x.size)

定义要拟合的函数：

def fit_msd2(params, x, data):
    A = params['A']
    B = params['B']
    c = params['c']

    model = A*x + B * (np.exp(-c*x)-1) 
    return model - data

params = Parameters()
params.add('A',   min=0, max = 100000)
params.add('B',   min=-100, max = 100000)
params.add('c',   min=0,max = 100000)

from scipy.optimize import dual_annealing
# do fit, here with the default leastsq algorithm

minner = Minimizer(fit_msd2, params, fcn_args=(x, y))
print(minner)
result = minner.minimize(method="dual_annealing")
print(result)
# calculate final result
final = x + result.residual
#print(final)

# write error report
report_fit(result)

fig, ax = plt.subplots(figsize = (10,10))
ax.grid()
ax.set_ylabel('$\Delta_{msd}$(t) [$\mu$m]', fontsize=18)
ax.set_xlabel('Time Lag $\Delta t$ [s]', fontsize=18)
ax.loglog(x, y, label = 'Synthetic data')
ax.loglog(x, final, 'r-', linewidth=3,label='Fit')
ax.legend(loc='best', fontsize = 12)

我认为配合不正常，我不明白是什么原因。我可以使用任何库或检查路由来评估它吗？ 如果没有，有人可以建议另一种方法来执行这种模拟退火拟合吗？

提前致谢

Answer 1

我最好的猜测是使用优化器的参数。由于您的损失似乎太大，优化器目前过早接受解决方案。由于lmfit 建立在scipy 之上，因此可以像这样传递优化器的参数：

# [...]
opt_args = {
  maxiter=5000,  # originally 1000
  initial_tempfloat=7500  # originally 5230
  accept=-2.  # originally -5.0
}
minner = Minimizer(fit_msd2, params, fcn_args=(x, y))
result = minner.minimize(method="dual_annealing", **opt_args)

上述参数化使优化器更具探索性，尤其是在优化开始时，它可以优化更长的时间。 dual_annealing 参数的完整列表可以在 here 找到。

这会带来改进吗？

另外，我建议使用平方你的错误：

def fit_msd2(params, x, data):
    A = params['A']
    B = params['B']
    c = params['c']

    model = A*x + B * (np.exp(-c*x)-1) 
    return (model - data)**2  # or np.pow(model - data, 2), depending on data type