Fit Two Dimensional Peaks

This example illustrates how to handle two-dimensional data with lmfit.

import matplotlib.pyplot as plt
import numpy as np
from scipy.interpolate import griddata

import lmfit
from lmfit.lineshapes import gaussian2d, lorentzian

Two-dimensional Gaussian

We start by considering a simple two-dimensional gaussian function, which depends on coordinates (x, y). The most general case of experimental data will be irregularly sampled and noisy. Let’s simulate some:

npoints = 10000
np.random.seed(2021)
x = np.random.rand(npoints)*10 - 4
y = np.random.rand(npoints)*5 - 3
z = gaussian2d(x, y, amplitude=30, centerx=2, centery=-.5, sigmax=.6, sigmay=.8)
z += 2*(np.random.rand(*z.shape)-.5)
error = np.sqrt(z+1)

To plot this, we can interpolate the data onto a grid.

X, Y = np.meshgrid(np.linspace(x.min(), x.max(), 100),
                   np.linspace(y.min(), y.max(), 100))
Z = griddata((x, y), z, (X, Y), method='linear', fill_value=0)

fig, ax = plt.subplots()
art = ax.pcolor(X, Y, Z, shading='auto')
plt.colorbar(art, ax=ax, label='z')
ax.set_xlabel('x')
ax.set_ylabel('y')
plt.show()

In this case, we can use a built-in model to fit

model = lmfit.models.Gaussian2dModel()
params = model.guess(z, x, y)
result = model.fit(z, x=x, y=y, params=params, weights=1/error)
lmfit.report_fit(result)

[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 87
    # data points      = 10000
    # variables        = 5
    chi-square         = 20618.1774
    reduced chi-square = 2.06284916
    Akaike info crit   = 7245.87992
    Bayesian info crit = 7281.93162
    R-squared          = 0.87508800
[[Variables]]
    amplitude:  27.4195833 +/- 0.65062974 (2.37%) (init = 16.51399)
    centerx:    1.99705425 +/- 0.01405864 (0.70%) (init = 1.940764)
    centery:   -0.49516158 +/- 0.01907800 (3.85%) (init = -0.5178641)
    sigmax:     0.54740777 +/- 0.01224965 (2.24%) (init = 1.666582)
    sigmay:     0.73300589 +/- 0.01617042 (2.21%) (init = 0.8332836)
    fwhmx:      1.28904676 +/- 0.02884573 (2.24%) == '2.3548200*sigmax'
    fwhmy:      1.72609693 +/- 0.03807842 (2.21%) == '2.3548200*sigmay'
    height:     10.8758308 +/- 0.35409534 (3.26%) == '0.1591549*amplitude/(max(1e-15, sigmax)*max(1e-15, sigmay))'
[[Correlations]] (unreported correlations are < 0.100)
    C(amplitude, sigmay) = +0.2464
    C(amplitude, sigmax) = +0.2314

To check the fit, we can evaluate the function on the same grid we used before and make plots of the data, the fit and the difference between the two.

fig, axs = plt.subplots(2, 2, figsize=(10, 10))

vmax = np.nanpercentile(Z, 99.9)

ax = axs[0, 0]
art = ax.pcolor(X, Y, Z, vmin=0, vmax=vmax, shading='auto')
plt.colorbar(art, ax=ax, label='z')
ax.set_title('Data')

ax = axs[0, 1]
fit = model.func(X, Y, **result.best_values)
art = ax.pcolor(X, Y, fit, vmin=0, vmax=vmax, shading='auto')
plt.colorbar(art, ax=ax, label='z')
ax.set_title('Fit')

ax = axs[1, 0]
fit = model.func(X, Y, **result.best_values)
art = ax.pcolor(X, Y, Z-fit, vmin=0, vmax=10, shading='auto')
plt.colorbar(art, ax=ax, label='z')
ax.set_title('Data - Fit')

for ax in axs.ravel():
    ax.set_xlabel('x')
    ax.set_ylabel('y')
axs[1, 1].remove()
plt.show()

Two-dimensional off-axis Lorentzian

We now go on to show a harder example, in which the peak has a Lorentzian profile and an off-axis anisotropic shape. This can be handled by applying a suitable coordinate transform and then using the lorentzian function that lmfit provides in the lineshapes module.

def lorentzian2d(x, y, amplitude=1., centerx=0., centery=0., sigmax=1., sigmay=1.,
                 rotation=0):
    """Return a two dimensional lorentzian.

    The maximum of the peak occurs at ``centerx`` and ``centery``
    with widths ``sigmax`` and ``sigmay`` in the x and y directions
    respectively. The peak can be rotated by choosing the value of ``rotation``
    in radians.
    """
    xp = (x - centerx)*np.cos(rotation) - (y - centery)*np.sin(rotation)
    yp = (x - centerx)*np.sin(rotation) + (y - centery)*np.cos(rotation)
    R = (xp/sigmax)**2 + (yp/sigmay)**2

    return 2*amplitude*lorentzian(R)/(np.pi*sigmax*sigmay)

Data can be simulated and plotted in the same way as we did before.

npoints = 10000
x = np.random.rand(npoints)*10 - 4
y = np.random.rand(npoints)*5 - 3
z = lorentzian2d(x, y, amplitude=30, centerx=2, centery=-.5, sigmax=.6,
                 sigmay=1.2, rotation=30*np.pi/180)
z += 2*(np.random.rand(*z.shape)-.5)
error = np.sqrt(z+1)

X, Y = np.meshgrid(np.linspace(x.min(), x.max(), 100),
                   np.linspace(y.min(), y.max(), 100))
Z = griddata((x, y), z, (X, Y), method='linear', fill_value=0)

fig, ax = plt.subplots()
ax.set_xlabel('x')
ax.set_ylabel('y')
art = ax.pcolor(X, Y, Z, shading='auto')
plt.colorbar(art, ax=ax, label='z')
plt.show()

To fit, create a model from the function. Don’t forget to tell lmfit that both x and y are independent variables. Keep in mind that lmfit will take the function keywords as default initial guesses in this case and that it will not know that certain parameters only make physical sense over restricted ranges. For example, peak widths should be positive and the rotation can be restricted over a quarter circle.

model = lmfit.Model(lorentzian2d, independent_vars=['x', 'y'])
params = model.make_params(amplitude=10, centerx=x[np.argmax(z)],
                           centery=y[np.argmax(z)])
params['rotation'].set(value=.1, min=0, max=np.pi/2)
params['sigmax'].set(value=1, min=0)
params['sigmay'].set(value=2, min=0)

result = model.fit(z, x=x, y=y, params=params, weights=1/error)
lmfit.report_fit(result)

[[Fit Statistics]]
    # fitting method   = leastsq
    # function evals   = 73
    # data points      = 10000
    # variables        = 6
    chi-square         = 11287.3823
    reduced chi-square = 1.12941588
    Akaike info crit   = 1223.00402
    Bayesian info crit = 1266.26606
    R-squared          = 0.85940121
[[Variables]]
    amplitude:  25.6417887 +/- 0.50569636 (1.97%) (init = 10)
    centerx:    2.00326033 +/- 0.01163397 (0.58%) (init = 2.051478)
    centery:   -0.49692376 +/- 0.01680902 (3.38%) (init = -0.478231)
    sigmax:     0.51074241 +/- 0.01027264 (2.01%) (init = 1)
    sigmay:     1.11741198 +/- 0.02190772 (1.96%) (init = 2)
    rotation:   0.48130689 +/- 0.01542118 (3.20%) (init = 0.1)
[[Correlations]] (unreported correlations are < 0.100)
    C(centerx, centery)  = +0.5767
    C(amplitude, sigmax) = +0.3342
    C(amplitude, sigmay) = +0.3046

The process of making plots to check it worked is the same as before.

fig, axs = plt.subplots(2, 2, figsize=(10, 10))

vmax = np.nanpercentile(Z, 99.9)

ax = axs[0, 0]
art = ax.pcolor(X, Y, Z, vmin=0, vmax=vmax, shading='auto')
plt.colorbar(art, ax=ax, label='z')
ax.set_title('Data')

ax = axs[0, 1]
fit = model.func(X, Y, **result.best_values)
art = ax.pcolor(X, Y, fit, vmin=0, vmax=vmax, shading='auto')
plt.colorbar(art, ax=ax, label='z')
ax.set_title('Fit')

ax = axs[1, 0]
fit = model.func(X, Y, **result.best_values)
art = ax.pcolor(X, Y, Z-fit, vmin=0, vmax=10, shading='auto')
plt.colorbar(art, ax=ax, label='z')
ax.set_title('Data - Fit')

for ax in axs.ravel():
    ax.set_xlabel('x')
    ax.set_ylabel('y')
axs[1, 1].remove()
plt.show()