Simulated annealing fit in Python

Question

I am trying to getting familiar with the non linear fitting procedure: dual-annealing. To do so I generated some synthetic data and try to fit over them a basic Furth formula, see the code below:

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 

generating synthetic data:

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)

Defining the function to fit:

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)

I believe the fit is not working properly and I cannot understand which are the reasons. Is there any library or check routing I can use in order to evaluate this? If not, can someone suggest another way to perform this simulated annealing fit?

Thanks in advance

Answer 1

My best guess would be to play around with the parameters of the optimizer. As it seems that your loss is too large, the optimizer currently accepts a solution too early. Since lmfit builds on scipy, it might be possible to pass the arguments of the optimizer as such:

# [...]
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)

The above parametrization makes the optimizers more explorative, especially in the beginning of the optimization and lets it optimize for longer. A full list of the dual_annealing parameters can be found here.

Does this lead to an improvement?

Also, I would suggest to use the square your error:

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