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
A1 2.03953400 0.30458039 (14.93%) 2 -inf inf True
A2 3.24628549 0.27392227 (8.44%) 3 -inf inf True
xc1 2.11179622 0.14262608 (6.75%) 2 -inf inf True
B 5.42073333 0.38660956 (7.13%) 5 -inf inf True
xc2 5.05045727 0.09355979 (1.85%) 5 -inf inf True
sigma1 0.82863194 0.18223396 (21.99%) 1 -inf inf True
sigma2 0.99878722 0.11821179 (11.84%) 1 -inf inf True
xw 3.51913583 0.43246203 (12.29%) 4 -inf inf True

Correlations (unreported correlations are < 0.100)

Bsigma1-0.6745
A2xw-0.6321
sigma2xw-0.5722
xc2sigma10.5537
A1A20.5293
A1xw-0.5152
xc1sigma2-0.4850
xc1xc20.4211
xc2sigma2-0.4149
sigma1sigma2-0.4054
A1B-0.3567
xc1B0.3213
A1sigma20.2476
Bxw-0.2475
Bxc2-0.2176
A1xc1-0.2079
A2sigma10.1766
Bsigma20.1662
A2B-0.1332
sigma1xw0.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-V3uHvC/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-V3uHvC/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
A1 2.16429205 0.17552276 (8.11%) 2 -inf inf True
xc1 2.13326031 0.14814919 (6.94%) 2 -inf inf True
B 5.17531388 0.29344865 (5.67%) 5 -inf inf True
xc2 5.09911007 0.07468792 (1.46%) 5 -inf inf True
sigma1 0.95831199 0.07211389 (7.53%) 1 -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
A1B-0.1919
A1sigma10.1428
xc1xw-0.1390


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

Gallery generated by Sphinx-Gallery