| Larch: Data Analysis Tools for X-ray Spectroscopy

# 13.7. Simplified Peak Fitting with `fit_peak()`

As shown in the previous sections, it is pretty simple to use Larch’s fitting mechanism to set up and perform fits to data. In fact, it is pretty commom to need to fit data to simple line-shapes, as when setting up an experiment. In an effort to make life easier, in addition to providing a fairly easy general-purpose fitting function, Larch provides a `fit_peak()` function to fit one dimensional data to one of several lineshapes.

_math.fit_peak(x, y, model, dy=None, background=None, step=None, negative=False, use_gamma=False)

  fit the data y at points x with a function described by model, optionally including uncertainties in y, a background polynomial, and a choice of step functions.

  Parameters:

    - x – array of values at which to calculate model

    - y – array of values for model to try to match

    - model – name of model for functional form to use. See Table of fit_peak models.

    - dy – optional array of values for uncertainty in y data to be matched.

    - background – name of background model to use. One of (case insensitive) None (no background), ‘constant’, ‘linear’, or ‘quadratic’. This is ignored when model is ‘linear’ or ‘quadratic’

    - step – name of step model to use for ‘step’ and ‘rectangle’ models. One of (case insensitive): ‘linear’, ‘erf’, or ‘atan’

    - negative – boolean for whether peak or steps are expected to go down (default `False`)

    - use_gamma – boolean for whether to use separate gamma parameter for ‘voigt’ model (default `False`).

  Returns:

    a Group containing fit parameters, best-fit and background functions, as detailed in Table of fit_peak group members.

> Table of `fit_peak()` models.

> > The following functional models are supported for fitting with `fit_peak()`.

> > model | Description | parameter names
> > --- | --- | ---
> > linear | \(y = mx + b\) | offset, slope
> > quadratic | \(y = a + bx + cx^2\) | offset, slope, quad
> > step | step function | height, center, width
> > rectangle | step up, step down | height, center1, width1, center2, width2
> > exponential | \(y = Ae^{-x/b}\) | amplitude, decay
> > gaussian | Gaussian | amplitude, center, sigma
> > lorentzian | Lorenztian | amplitude, center, sigma
> > voigt | Voigt (with gamma=sigma) | amplitude, center, sigma

The sigma and width parameters will have a minimum value of 0.

For all models except linear and quadratic, an optional background can be included, which (depending on the form chosen) will add parameters named bkg_offset, bkg_slope, and bkg_quad.

> Table of members of group returned by `fit_peak()`.

> > member | Description
> > --- | ---
> > x | copied from input x data
> > y | copied from input y data
> > model | name of model used
> > background | name of background models used
> > step | name of step model used
> > params | Group of fit parameters, as sent to and from `minimize()`
> > fit | final best fit to y data, including background if used.
> > fit_init | initial guess of fit to y data, including background if used.
> > bkg | final background function, if used.
> > bkg_init | initial background function, if used.

## 13.7.1. Example: Fitting a Gaussian + background with `fit_peak()`

As in the Example in the previous section, we make a simple mock data set and fit a Gaussian function to it. Here we also add a linear background, and do the whole fit with a single function, instead of a dozen or so lines of code used before:

    ## examples/fitting/doc_example_fitpeak1.lar
    # create mock data
    x = linspace(0, 100, 51)
    noise = random.normal(size=len(x), scale=0.1)
    y = 8.0 - x*0.025 + noise
    y = y + voigt(x, amplitude=89.0, center=44.25, sigma=4.75, gamma=2.0)

    # fit to Gaussian with linear background
    myfit = fit_peak(x, y, 'gaussian', background='linear')

    plot(myfit.x, myfit.y, marker='+', label='data',
         xlabel='x', ylabel='y', show_legend=True, new=True)
    plot(myfit.x, myfit.fit_init, label='init')
    plot(myfit.x, myfit.fit, label='best fit')

    print( fit_report(myfit, min_correl=0.3))

    ## end examples/fitting/doc_example_fitpeak1.lar

Here we first create mock data that’s fairly noisy, and ask it to be fit to Gaussian, including a linear background, with a single command. We then plot the results and print the report of parameters. Note that the fit is pretty good at finding the peak center and shape, even though the model is not strictly correct.

The printed output from `fit_report(myfit)` will look like this:

    [[Model]]
        (Model(gaussian) + Model(linear, prefix='bkg_'))
    [[Fit Statistics]]
        # function evals   = 33
        # data points      = 51
        # variables        = 5
        chi-square         = 0.62209
        reduced chi-square = 0.013524
        Akaike info crit   = -214.73
        Bayesian info crit = -205.07
    [[Variables]]
        amplitude:       75.4295207 +/- 1.085962 (1.44%) (init= 143.506)
        bkg_intercept:   8.10849304 +/- 0.035243 (0.43%) (init= 0)
        bkg_slope:      -0.02506582 +/- 0.000562 (2.24%) (init= 0)
        center:          44.3598293 +/- 0.079202 (0.18%) (init= 43)
        fwhm:            13.3776066 +/- 0.198318 (1.48%)  == '2.3548200*sigma'
        height:          5.29700921 +/- 0.065021 (1.23%)  == 0.3989423*amplitude/max(1.e-15, sigma)'
        sigma:           5.68094659 +/- 0.084218 (1.48%) (init= 7)
    [[Correlations]] (unreported correlations are <  0.300)
        C(bkg_intercept, bkg_slope)  = -0.835
        C(amplitude, sigma)          =  0.647
        C(amplitude, bkg_intercept)  = -0.402

And the plot of data and fit will look like this:

Figure 13.7.1.1 Simple fit to mock data using a Gaussian model and a linear background with the `fit_peak()` function.

Although the fit is quite good, the model is probably imperfect, and using a Voigt function to fit to this data would give better results. The main point here is not just that the fit is good, but that it was accomplished with a single line of code.

## 13.7.2. Example: Fitting a Rectangular function with `fit_peak()`

Here, we simulate data that might represent a line-up scan at a beamline. The fit is to a function that takes a step up and a step down.

    ## examples/fitting/doc_example_fitpeak2.lar
    # create mock rectangular data
    def make_rect_data(x, noise_scale=0.01):
        noise = random.normal(size=len(x), scale=noise_scale)
        y = 6.0 - x*0.001 + noise
        x0 = 30
        dx = 7
        sig = dx*0.6
        y = y + gaussian(x, amplitude=68, center=x0+0*dx, sigma=1.2*sig)
        y = y + gaussian(x, amplitude=60, center=x0+1*dx, sigma=1.1*sig)
        y = y + gaussian(x, amplitude=71, center=x0+2*dx, sigma=sig)
        y = y + gaussian(x, amplitude=75, center=x0+3*dx, sigma=sig)
        y = y + gaussian(x, amplitude=70, center=x0+4*dx, sigma=sig)
        y = y + gaussian(x, amplitude=65, center=x0+5*dx, sigma=sig)
        y = y + gaussian(x, amplitude=72, center=x0+6*dx, sigma=sig)
        y = y + gaussian(x, amplitude=61, center=x0+7*dx, sigma=0.9*sig)
        return y
    #enddef

    x = linspace(0, 160, 321)
    y = make_rect_data(x, noise_scale=0.2)

    # fit to rectangle shape
    myfit = fit_peak(x, y, 'rectangle', background='constant', step='erf')

    plot(myfit.x, myfit.y, marker='+', label='data',
         xlabel='x', ylabel='y', show_legend=True, new=True)
    plot(myfit.x, myfit.fit, label='best fit')

    print( fit_report(myfit, min_correl=0.3))

    ## end examples/fitting/doc_example_fitpeak2.lar

Again, we first create mock data that’s fairly noisy. The data is clearly not exactly rectangular, but we ask it to be fit to a rectangular function plus a linear background. The printed output from `fit_report(myfit)` will look like this:

    [[Model]]
        (Model(rectangle, form='erf') + Model(constant, prefix='bkg_'))
    [[Fit Statistics]]
        # function evals   = 67
        # data points      = 321
        # variables        = 6
        chi-square         = 23.405
        reduced chi-square = 0.074301
        Akaike info crit   = -828.54
        Bayesian info crit = -805.91
    [[Variables]]
        amplitude:   9.95178246 +/- 0.038770 (0.39%) (init= 11.34094)
        bkg_c:       5.90041256 +/- 0.020680 (0.35%) (init= 0)
        center1:     27.1354829 +/- 0.095553 (0.35%) (init= 26.5)
        center2:     82.0689673 +/- 0.071565 (0.09%) (init= 81.5)
        midpoint:    54.6022251 +/- 0.055873 (0.10%)  == '(center1+center2)/2.0'
        sigma1:      8.14570023 +/- 0.186184 (2.29%) (init= 32)
        sigma2:      4.88414083 +/- 0.140778 (2.88%) (init= 32)
    [[Correlations]] (unreported correlations are <  0.300)
        C(amplitude, bkg_c)          = -0.566
        C(amplitude, sigma1)         =  0.341

Figure 13.7.2.1 Simple fit to mock data to a rectangular function and a linear background using the `fit_peak()` function.

Again, the principle point here is not how well the rectangular model matches the actual data here, but how simply one can model data to a selection of simple shapes. Such fits can be very useful for preliminary data visualization and analysis. Of course, one should be cautious about treating the results of such an automated approach as a final and best analysis of any data.

## 13.7.3. Using lmfit.Model

Note that the `fit_peak()` function gives a simple wrapping of lmfit.Model. For a wider selection of builtin Models and more sophisticated model building including adding bounds and constraints between parameters one can import and use lmfit.Model directly with larch. As an example, the above fit can be replicated with:

    larch> from lmfit.models import RectangleModel, ConstantModel
    larch> model = RectangleModel(form='erf') + ConstantModel(prefix='bkg_')
    larch> params = model.make_params(bkg_c=0, sigma1=10, sigma2=10,
                                      center1=20, center20=80, amplitude=10)
    larch> out = model.fit(y, params, x=x)
    larch> print(out.fit_report(sorted_pars=True))

Here, a model is built as the sum of two components: a rectangle and a constant background. lmfit.Parameters are made from this composite model and the parameter names, giving initial values for each parameter. The model can then be used to fit the data `y` with the parameters `params`, and the independent variable `x`. While directly using lmfit.Model does require providing initial values for all parameters, it also gives complete access to the parameters, and allows building more complex models.