- Contents
- Examples
- Installation
- FAQ
- Support
- Develop

Note

Click here to download the full example code

# 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')

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

Model(_lambdifygenerated)

## Fit Statistics

fitting method | leastsq |
# function evals | 64 |
# data points | 40 |
# variables | 8 |
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)

B | sigma1 | -0.6745
A2 | xw | -0.6321
sigma2 | xw | -0.5722
xc2 | sigma1 | 0.5537
A1 | A2 | 0.5293
A1 | xw | -0.5152
xc1 | sigma2 | -0.4850
xc1 | xc2 | 0.4211
xc2 | sigma2 | -0.4149
sigma1 | sigma2 | -0.4054
A1 | B | -0.3567
xc1 | B | 0.3213
A1 | sigma2 | 0.2476
B | xw | -0.2475
B | xc2 | -0.2176
A1 | xc1 | -0.2079
A2 | sigma1 | 0.1766
B | sigma2 | 0.1662
A2 | B | -0.1332
sigma1 | xw | 0.1048
A2 | xc2 | 0.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

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 method | leastsq |
# function evals | 29 |
# data points | 40 |
# variables | 6 |
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)

A1 | xw | -0.6832
xc1 | B | 0.6816
B | sigma1 | -0.4525
xc1 | sigma1 | -0.4283
xc1 | xc2 | 0.4184
sigma1 | xw | -0.4158
B | xc2 | 0.3071
xc2 | xw | -0.2282
B | xw | -0.1931
A1 | B | -0.1919
A1 | sigma1 | 0.1428
xc1 | xw | -0.1390

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

Gallery generated by Sphinx-Gallery