Fitting lmfit.Models to xarray objects

xarray-lmfit adds some methods to xarray objects that allows you to fit data with lmfit models: xarray.DataArray.xlm.modelfit() and xarray.Dataset.xlm.modelfit(), depending on whether you want to fit a single DataArray or multiple DataArrays in a Dataset.

Hint

The syntax of the accessors are similar to the xarray native methods xarray.DataArray.curvefit() and xarray.Dataset.curvefit().

First, let us generate a Gaussian peak on a linear background.

import numpy as np
import xarray as xr
import matplotlib.pyplot as plt
import lmfit

# Generate Gaussian peak on linear background
x = np.linspace(0, 10, 50)
y = -0.1 * x + 2 + 3 * np.exp(-((x - 5) ** 2) / (2 * 1**2))

# Add some noise with fixed seed for reproducibility
rng = np.random.default_rng(5)
yerr = np.full_like(x, 0.3)
y = rng.normal(y, yerr)

y_arr = xr.DataArray(y, dims=("x",), coords={"x": x})

y_arr

<xarray.DataArray (x: 50)> Size: 400B
array([1.75943175, 1.58231452, 1.88475436, 2.06510701, 2.25965595,
       1.93196286, 1.71416562, 1.62680659, 2.07170706, 2.32691224,
       1.91538585, 1.47294052, 1.58349782, 2.40570597, 2.07715715,
       1.6322706 , 2.31469126, 2.23390379, 2.68841232, 3.06370507,
       3.34563648, 4.0619544 , 4.16597499, 4.21645875, 4.61751752,
       4.7231746 , 3.83915663, 4.00584093, 3.45885951, 3.32401046,
       2.59879139, 2.61809975, 2.26538341, 1.90156002, 1.449778  ,
       1.46886389, 1.12891382, 0.95421068, 1.35935443, 0.9089135 ,
       1.5549131 , 1.38859024, 0.54860997, 1.20649644, 0.77261759,
       1.09202446, 1.07451202, 0.44436656, 0.95041178, 0.92327417])
Coordinates:
  * x        (x) float64 400B 0.0 0.2041 0.4082 0.6122 ... 9.592 9.796 10.0

xarray.DataArray

• x: 50

• 1.759 1.582 1.885 2.065 2.26 ... 1.092 1.075 0.4444 0.9504 0.9233

  array([1.75943175, 1.58231452, 1.88475436, 2.06510701, 2.25965595,
         1.93196286, 1.71416562, 1.62680659, 2.07170706, 2.32691224,
         1.91538585, 1.47294052, 1.58349782, 2.40570597, 2.07715715,
         1.6322706 , 2.31469126, 2.23390379, 2.68841232, 3.06370507,
         3.34563648, 4.0619544 , 4.16597499, 4.21645875, 4.61751752,
         4.7231746 , 3.83915663, 4.00584093, 3.45885951, 3.32401046,
         2.59879139, 2.61809975, 2.26538341, 1.90156002, 1.449778  ,
         1.46886389, 1.12891382, 0.95421068, 1.35935443, 0.9089135 ,
         1.5549131 , 1.38859024, 0.54860997, 1.20649644, 0.77261759,
         1.09202446, 1.07451202, 0.44436656, 0.95041178, 0.92327417])

• Coordinates: (1)
  • x
    (x)
    float64
    0.0 0.2041 0.4082 ... 9.796 10.0

    array([ 0.      ,  0.204082,  0.408163,  0.612245,  0.816327,  1.020408,
            1.22449 ,  1.428571,  1.632653,  1.836735,  2.040816,  2.244898,
            2.44898 ,  2.653061,  2.857143,  3.061224,  3.265306,  3.469388,
            3.673469,  3.877551,  4.081633,  4.285714,  4.489796,  4.693878,
            4.897959,  5.102041,  5.306122,  5.510204,  5.714286,  5.918367,
            6.122449,  6.326531,  6.530612,  6.734694,  6.938776,  7.142857,
            7.346939,  7.55102 ,  7.755102,  7.959184,  8.163265,  8.367347,
            8.571429,  8.77551 ,  8.979592,  9.183673,  9.387755,  9.591837,
            9.795918, 10.      ])

• Indexes: (1)
  • x
    PandasIndex

    PandasIndex(Index([                0.0, 0.20408163265306123, 0.40816326530612246,
            0.6122448979591837,  0.8163265306122449,  1.0204081632653061,
            1.2244897959183674,  1.4285714285714286,  1.6326530612244898,
             1.836734693877551,  2.0408163265306123,  2.2448979591836737,
            2.4489795918367347,  2.6530612244897958,   2.857142857142857,
            3.0612244897959187,  3.2653061224489797,  3.4693877551020407,
             3.673469387755102,  3.8775510204081636,   4.081632653061225,
             4.285714285714286,  4.4897959183673475,  4.6938775510204085,
            4.8979591836734695,  5.1020408163265305,  5.3061224489795915,
             5.510204081632653,   5.714285714285714,   5.918367346938775,
             6.122448979591837,   6.326530612244898,   6.530612244897959,
              6.73469387755102,   6.938775510204081,   7.142857142857143,
             7.346938775510204,   7.551020408163265,   7.755102040816327,
             7.959183673469388,    8.16326530612245,    8.36734693877551,
             8.571428571428571,   8.775510204081632,   8.979591836734695,
             9.183673469387756,   9.387755102040817,   9.591836734693878,
             9.795918367346939,                10.0],
          dtype='float64', name='x'))

• Attributes: (0)

Here, y_arr is the DataArray that contains both the values we want to fit, along with the independent variable x as a coordinate.

After importing xarray_lmfit, we can use xarray.DataArray.xlm.modelfit() to fit the data to a model.

import xarray_lmfit

model = lmfit.models.GaussianModel() + lmfit.models.LinearModel()

params = model.make_params(slope=-0.1, center=5.0, sigma={"value": 0.1, "min": 0})

result = y_arr.xlm.modelfit("x", model=model, params=params)

The fit result Dataset

xarray.DataArray.xlm.modelfit() returns a xarray.Dataset including the best-fit parameters and the fit statistics:

result

<xarray.Dataset> Size: 2kB
Dimensions:                (param: 5, cov_i: 5, cov_j: 5, fit_stat: 9, x: 50)
Coordinates:
  * x                      (x) float64 400B 0.0 0.2041 0.4082 ... 9.796 10.0
  * param                  (param) <U9 180B 'amplitude' 'center' ... 'intercept'
  * fit_stat               (fit_stat) <U8 288B 'nfev' 'nvarys' ... 'rsquared'
  * cov_i                  (cov_i) <U9 180B 'amplitude' 'center' ... 'intercept'
  * cov_j                  (cov_j) <U9 180B 'amplitude' 'center' ... 'intercept'
Data variables:
    modelfit_results       object 8B <lmfit.model.ModelResult object at 0x77d...
    modelfit_coefficients  (param) float64 40B 7.166 4.992 0.943 -0.1138 2.006
    modelfit_stderr        (param) float64 40B 0.3762 0.04221 ... 0.08441
    modelfit_covariance    (cov_i, cov_j) float64 200B 0.1415 ... 0.007126
    modelfit_stats         (fit_stat) float64 72B 43.0 5.0 ... -118.4 0.9468
    modelfit_data          (x) float64 400B 1.759 1.582 1.885 ... 0.9504 0.9233
    modelfit_best_fit      (x) float64 400B 2.006 1.983 1.96 ... 0.8915 0.8683

xarray.Dataset

• Dimensions:
  • param: 5
  • cov_i: 5
  • cov_j: 5
  • fit_stat: 9
  • x: 50
• Coordinates: (5)
  • x
    (x)
    float64
    0.0 0.2041 0.4082 ... 9.796 10.0

    array([ 0.      ,  0.204082,  0.408163,  0.612245,  0.816327,  1.020408,
            1.22449 ,  1.428571,  1.632653,  1.836735,  2.040816,  2.244898,
            2.44898 ,  2.653061,  2.857143,  3.061224,  3.265306,  3.469388,
            3.673469,  3.877551,  4.081633,  4.285714,  4.489796,  4.693878,
            4.897959,  5.102041,  5.306122,  5.510204,  5.714286,  5.918367,
            6.122449,  6.326531,  6.530612,  6.734694,  6.938776,  7.142857,
            7.346939,  7.55102 ,  7.755102,  7.959184,  8.163265,  8.367347,
            8.571429,  8.77551 ,  8.979592,  9.183673,  9.387755,  9.591837,
            9.795918, 10.      ])

  • param
    (param)
    <U9
    'amplitude' ... 'intercept'

    array(['amplitude', 'center', 'sigma', 'slope', 'intercept'], dtype='<U9')

  • fit_stat
    (fit_stat)
    <U8
    'nfev' 'nvarys' ... 'rsquared'

    array(['nfev', 'nvarys', 'ndata', 'nfree', 'chisqr', 'redchi', 'aic', 'bic',
           'rsquared'], dtype='<U8')

  • cov_i
    (cov_i)
    <U9
    'amplitude' ... 'intercept'

    array(['amplitude', 'center', 'sigma', 'slope', 'intercept'], dtype='<U9')

  • cov_j
    (cov_j)
    <U9
    'amplitude' ... 'intercept'

    array(['amplitude', 'center', 'sigma', 'slope', 'intercept'], dtype='<U9')

• Data variables: (7)
  • modelfit_results
    ()
    object
    <lmfit.model.ModelResult object ...

    array(<lmfit.model.ModelResult object at 0x77db575306b0>, dtype=object)

  • modelfit_coefficients
    (param)
    float64
    7.166 4.992 0.943 -0.1138 2.006

    array([ 7.16590435,  4.99152727,  0.9430048 , -0.11380431,  2.0063381 ])

  • modelfit_stderr
    (param)
    float64
    0.3762 0.04221 ... 0.01319 0.08441

    array([0.37616176, 0.04220625, 0.04687009, 0.01318512, 0.08441486])

  • modelfit_covariance
    (cov_i, cov_j)
    float64
    0.1415 -1.2e-05 ... 0.007126

    array([[ 1.41497670e-01, -1.20037799e-05,  1.24134218e-02,
             1.44689850e-05, -1.39390991e-02],
           [-1.20037799e-05,  1.78136789e-03, -1.10018585e-06,
            -1.44229983e-04,  7.22326264e-04],
           [ 1.24134218e-02, -1.10018585e-06,  2.19680556e-03,
             1.27311822e-06, -1.22287791e-03],
           [ 1.44689850e-05, -1.44229983e-04,  1.27311822e-06,
             1.73847451e-04, -8.70655213e-04],
           [-1.39390991e-02,  7.22326264e-04, -1.22287791e-03,
            -8.70655213e-04,  7.12586920e-03]])

  • modelfit_stats
    (fit_stat)
    float64
    43.0 5.0 50.0 ... -118.4 0.9468

    array([ 4.30000000e+01,  5.00000000e+00,  5.00000000e+01,  4.50000000e+01,
            3.16476712e+00,  7.03281583e-02, -1.27997177e+02, -1.18437062e+02,
            9.46775170e-01])

  • modelfit_data
    (x)
    float64
    1.759 1.582 1.885 ... 0.9504 0.9233

    array([1.75943175, 1.58231452, 1.88475436, 2.06510701, 2.25965595,
           1.93196286, 1.71416562, 1.62680659, 2.07170706, 2.32691224,
           1.91538585, 1.47294052, 1.58349782, 2.40570597, 2.07715715,
           1.6322706 , 2.31469126, 2.23390379, 2.68841232, 3.06370507,
           3.34563648, 4.0619544 , 4.16597499, 4.21645875, 4.61751752,
           4.7231746 , 3.83915663, 4.00584093, 3.45885951, 3.32401046,
           2.59879139, 2.61809975, 2.26538341, 1.90156002, 1.449778  ,
           1.46886389, 1.12891382, 0.95421068, 1.35935443, 0.9089135 ,
           1.5549131 , 1.38859024, 0.54860997, 1.20649644, 0.77261759,
           1.09202446, 1.07451202, 0.44436656, 0.95041178, 0.92327417])

  • modelfit_best_fit
    (x)
    float64
    2.006 1.983 1.96 ... 0.8915 0.8683

    array([2.0063406 , 1.98312041, 1.95990985, 1.93672489, 1.91360446,
           1.89063868, 1.86802457, 1.84616915, 1.82586501, 1.80856415,
           1.79676144, 1.79446144, 1.80763442, 1.8444757 , 1.91519567,
           2.03103699, 2.20230617, 2.43546092, 2.72969479, 3.07389277,
           3.44510291, 3.80957292, 4.12683745, 4.3564062 , 4.46560962,
           4.43652413, 4.26995411, 3.98525072, 3.6160237 , 3.20305677,
           2.78648377, 2.39927395, 2.06339527, 1.78901996, 1.57622315,
           1.41809979, 1.30416605, 1.22320425, 1.16514639, 1.12198239,
           1.08791937, 1.05909659, 1.03312278, 1.00861383, 0.98481809,
           0.96135149, 0.93802899, 0.91476638, 0.89152741, 0.86829731])

• Indexes: (5)
  • x
    PandasIndex

    PandasIndex(Index([                0.0, 0.20408163265306123, 0.40816326530612246,
            0.6122448979591837,  0.8163265306122449,  1.0204081632653061,
            1.2244897959183674,  1.4285714285714286,  1.6326530612244898,
             1.836734693877551,  2.0408163265306123,  2.2448979591836737,
            2.4489795918367347,  2.6530612244897958,   2.857142857142857,
            3.0612244897959187,  3.2653061224489797,  3.4693877551020407,
             3.673469387755102,  3.8775510204081636,   4.081632653061225,
             4.285714285714286,  4.4897959183673475,  4.6938775510204085,
            4.8979591836734695,  5.1020408163265305,  5.3061224489795915,
             5.510204081632653,   5.714285714285714,   5.918367346938775,
             6.122448979591837,   6.326530612244898,   6.530612244897959,
              6.73469387755102,   6.938775510204081,   7.142857142857143,
             7.346938775510204,   7.551020408163265,   7.755102040816327,
             7.959183673469388,    8.16326530612245,    8.36734693877551,
             8.571428571428571,   8.775510204081632,   8.979591836734695,
             9.183673469387756,   9.387755102040817,   9.591836734693878,
             9.795918367346939,                10.0],
          dtype='float64', name='x'))

  • param
    PandasIndex

    PandasIndex(Index(['amplitude', 'center', 'sigma', 'slope', 'intercept'], dtype='object', name='param'))

  • fit_stat
    PandasIndex

    PandasIndex(Index(['nfev', 'nvarys', 'ndata', 'nfree', 'chisqr', 'redchi', 'aic', 'bic',
           'rsquared'],
          dtype='object', name='fit_stat'))

  • cov_i
    PandasIndex

    PandasIndex(Index(['amplitude', 'center', 'sigma', 'slope', 'intercept'], dtype='object', name='cov_i'))

  • cov_j
    PandasIndex

    PandasIndex(Index(['amplitude', 'center', 'sigma', 'slope', 'intercept'], dtype='object', name='cov_j'))

• Attributes: (0)

Let’s take a closer look at the data variables in the resulting Dataset.

• modelfit_results contains the underlying lmfit.model.ModelResult object from the fit.

• modelfit_coefficients and modelfit_stderr contain the best-fit coefficients and their errors, respectively.

• modelfit_stats contains the goodness-of-fit statistics.

When called on a Dataset instead of a DataArray, these variables will be prefixed with the name of the data variable they correspond to.

It may not be immediately obvious how this is useful, but the true power of the accessor comes from its ability to utilize xarray’s powerful broadcasting capabilities, as described in the next section.