Fit Using differential_evolution Algorithm

This example compares the “leastsq” and “differential_evolution” algorithms on a fairly simple problem.

example diffev

Out:

# Fit using leastsq:
[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 65
    # data points      = 101
    # variables        = 4
    chi-square         = 21.7961792
    reduced chi-square = 0.22470288
    Akaike info crit   = -146.871969
    Bayesian info crit = -136.411487
[[Variables]]
    offset:  0.96333090 +/- 0.04735921 (4.92%) (init = 2)
    omega:   3.98700821 +/- 0.02079710 (0.52%) (init = 3.3)
    amp:     1.80253568 +/- 0.19401987 (10.76%) (init = 2.5)
    decay:   5.76279886 +/- 1.04073317 (18.06%) (init = 1)
[[Correlations]] (unreported correlations are < 0.100)
    C(amp, decay) = -0.755


# Fit using differential_evolution:
[[Fit Statistics]]
    # fitting method   = differential_evolution
    # function evals   = 1425
    # data points      = 101
    # variables        = 4
    chi-square         = 21.7961792
    reduced chi-square = 0.22470288
    Akaike info crit   = -146.871969
    Bayesian info crit = -136.411487
##  Warning: uncertainties could not be estimated:
    this fitting method does not natively calculate uncertainties
    and numdifftools is not installed for lmfit to do this. Use
    `pip install numdifftools` for lmfit to estimate uncertainties
    with this fitting method.
[[Variables]]
    offset:  0.96333133 (init = 2)
    omega:   3.98700854 (init = 3.3)
    amp:     1.80252613 (init = 2.5)
    decay:   5.76284516 (init = 1)

import matplotlib.pyplot as plt
import numpy as np

import lmfit

np.random.seed(2)
x = np.linspace(0, 10, 101)

# Setup example
decay = 5
offset = 1.0
amp = 2.0
omega = 4.0

y = offset + amp*np.sin(omega*x) * np.exp(-x/decay)
yn = y + np.random.normal(size=y.size, scale=0.450)


def resid(params, x, ydata):
    decay = params['decay'].value
    offset = params['offset'].value
    omega = params['omega'].value
    amp = params['amp'].value

    y_model = offset + amp * np.sin(x*omega) * np.exp(-x/decay)
    return y_model - ydata


params = lmfit.Parameters()
params.add('offset', 2.0, min=0, max=10.0)
params.add('omega', 3.3, min=0, max=10.0)
params.add('amp', 2.5, min=0, max=10.0)
params.add('decay', 1.0, min=0, max=10.0)

o1 = lmfit.minimize(resid, params, args=(x, yn), method='leastsq')
print("# Fit using leastsq:")
lmfit.report_fit(o1)

o2 = lmfit.minimize(resid, params, args=(x, yn), method='differential_evolution')
print("\n\n# Fit using differential_evolution:")
lmfit.report_fit(o2)

plt.plot(x, yn, 'ko', lw=2)
plt.plot(x, yn+o1.residual, 'r-', lw=2)
plt.plot(x, yn+o2.residual, 'b--', lw=2)
plt.legend(['data', 'leastsq', 'diffev'], loc='upper left')
plt.show()

Total running time of the script: ( 0 minutes 0.413 seconds)

Gallery generated by Sphinx-Gallery