Building a lmfit model with SymPy

SymPy is a Python library for symbolic mathematics. It can be very useful to build a model with sympy and then use that apply that model to the data with lmfit. This example shows how to do that. Notice, that this example requires both sympy and matplotlib.

import matplotlib.pyplot as plt
import numpy as np
import sympy
from sympy.parsing import sympy_parser

import lmfit

np.random.seed(1)

Instead of creating the sympy-sybols explicitly and building an expression with them, we will use the sympy parser.

gauss_peak1 = sympy_parser.parse_expr('A1*exp(-(x-xc1)**2/(2*sigma1**2))')
gauss_peak2 = sympy_parser.parse_expr('A2*exp(-(x-xc2)**2/(2*sigma2**2))')
exp_back = sympy_parser.parse_expr('B*exp(-x/xw)')

model_list = sympy.Array((gauss_peak1, gauss_peak2, exp_back))
model = sum(model_list)
model

Out:

A1*exp(-(x - xc1)**2/(2*sigma1**2)) + A2*exp(-(x - xc2)**2/(2*sigma2**2)) + B*exp(-x/xw)

We are using sympys lambdify function to make a function from the model expressions. We use these functions to generate some fake data.

model_list_func = sympy.lambdify(list(model_list.free_symbols), model_list)
model_func = sympy.lambdify(list(model.free_symbols), model)

x = np.linspace(0, 10, 40)
param_values = dict(x=x, A1=2, sigma1=1, sigma2=1, A2=3,
                    xc1=2, xc2=5, xw=4, B=5)
y = model_func(**param_values)
yi = model_list_func(**param_values)
yn = y + np.random.randn(y.size)*0.4

plt.plot(x, yn, 'o', zorder=1.9, ms=3)
plt.plot(x, y, lw=3)
for c in yi:
    plt.plot(x, c, lw=1, c='0.7')
example sympy

Next, we will just create a lmfit model from the function and fit the data.

lm_mod = lmfit.Model(model_func, independent_vars=('x'))
res = lm_mod.fit(data=yn, **param_values)
res.plot_fit()
plt.plot(x, y, label='true')
plt.legend()

res
Model(_lambdifygenerated)

Model

Model(_lambdifygenerated)

Fit Statistics

fitting methodleastsq
# function evals64
# data points40
# variables8
chi-square 5.39675674
reduced chi-square 0.16864865
Akaike info crit.-64.1232514
Bayesian info crit.-50.6122157

Variables

name value standard error relative error initial value min max vary
xc1 2.11179622 0.14262608 (6.75%) 2 -inf inf True
sigma2 0.99878722 0.11821179 (11.84%) 1 -inf inf True
A2 3.24628549 0.27392227 (8.44%) 3 -inf inf True
B 5.42073333 0.38660956 (7.13%) 5 -inf inf True
xw 3.51913583 0.43246203 (12.29%) 4 -inf inf True
xc2 5.05045727 0.09355979 (1.85%) 5 -inf inf True
A1 2.03953400 0.30458039 (14.93%) 2 -inf inf True
sigma1 0.82863194 0.18223396 (21.99%) 1 -inf inf True

Correlations (unreported correlations are < 0.100)

Bsigma1-0.6745
A2xw-0.6321
sigma2xw-0.5722
xc2sigma10.5537
A2A10.5293
xwA1-0.5152
xc1sigma2-0.4850
xc1xc20.4211
sigma2xc2-0.4149
sigma2sigma1-0.4054
BA1-0.3567
xc1B0.3213
sigma2A10.2476
Bxw-0.2475
Bxc2-0.2176
xc1A1-0.2079
A2sigma10.1766
sigma2B0.1662
A2B-0.1332
xwsigma10.1048
A2xc20.1032


The nice thing of using sympy is that we can easily modify our fit function. Let’s assume we know that the width of both gaussians is identical. Simliary, we assume that the ratio between both gaussians is fixed to 3:2 for some reason. Both can be expressed by just substituting the variables.

model2 = model.subs('sigma2', 'sigma1').subs('A2', '3/2*A1')
model2_func = sympy.lambdify(list(model2.free_symbols), model2)
lm_mod = lmfit.Model(model2_func, independent_vars=('x'))
param2_values = dict(x=x, A1=2, sigma1=1, A2=3, xc1=2, xc2=5, xw=4, B=5)
res2 = lm_mod.fit(data=yn, **param_values)
res2.plot_fit()
plt.plot(x, y, label='true')
plt.legend()

res2
Model(_lambdifygenerated)

Out:

/build/lmfit-py-4Y2iDh/lmfit-py-1.0.2/lmfit/model.py:968: UserWarning: The keyword argument sigma2 does not match any arguments of the model function. It will be ignored.
  warnings.warn("The keyword argument %s does not " % name +
/build/lmfit-py-4Y2iDh/lmfit-py-1.0.2/lmfit/model.py:968: UserWarning: The keyword argument A2 does not match any arguments of the model function. It will be ignored.
  warnings.warn("The keyword argument %s does not " % name +

Model

Model(_lambdifygenerated)

Fit Statistics

fitting methodleastsq
# function evals29
# data points40
# variables6
chi-square 5.49246202
reduced chi-square 0.16154300
Akaike info crit.-67.4201137
Bayesian info crit.-57.2868370

Variables

name value standard error relative error initial value min max vary
xc1 2.13326031 0.14814919 (6.94%) 2 -inf inf True
B 5.17531388 0.29344865 (5.67%) 5 -inf inf True
sigma1 0.95831199 0.07211389 (7.53%) 1 -inf inf True
xc2 5.09911007 0.07468792 (1.46%) 5 -inf inf True
A1 2.16429205 0.17552276 (8.11%) 2 -inf inf True
xw 3.56835239 0.40042300 (11.22%) 4 -inf inf True

Correlations (unreported correlations are < 0.100)

A1xw-0.6832
xc1B0.6816
Bsigma1-0.4525
xc1sigma1-0.4283
xc1xc20.4184
sigma1xw-0.4158
Bxc20.3071
xc2xw-0.2282
Bxw-0.1931
BA1-0.1919
sigma1A10.1428
xc1xw-0.1390


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

Gallery generated by Sphinx-Gallery