LM_FIT
Overview
The LM_FIT function performs nonlinear least-squares curve fitting using the lmfit library’s extensive collection of built-in models. It enables users to fit experimental data to mathematical models ranging from simple polynomials to complex peak shapes, with support for composite models that combine multiple functions.
The function leverages the Levenberg-Marquardt algorithm (LMA) by default, a widely-used optimization method for solving nonlinear least squares problems. LMA interpolates between the Gauss-Newton algorithm and gradient descent, offering robust convergence even when initial parameter estimates are far from optimal. The algorithm minimizes the sum of squared residuals:
S(\beta) = \sum_{i=1}^{m} [y_i - f(x_i, \beta)]^2where y_i are the observed values, f(x_i, \beta) is the model function, and \beta represents the parameters to be optimized.
The available models span several categories: peak-like models (Gaussian, Lorentzian, Voigt, Pseudo-Voigt, and various skewed/asymmetric distributions), polynomial models (linear, quadratic, and higher-order polynomials), exponential and power-law models, periodic models (sine waves, damped oscillators), and step-like models (step functions, Fermi distributions). Each model includes automatic parameter guessing via the guess() method and sensible default constraints.
A key feature is composite model fitting, which allows combining two models to fit complex datasets. For example, a Gaussian peak can be added to a linear baseline, or two overlapping peaks can be modeled simultaneously. When models are combined, parameters are prefixed (e.g., m1_ and m2_) to avoid naming conflicts.
The function returns fitted parameter values along with their standard errors, enabling uncertainty quantification in the fitted results. Alternative fitting methods including Trust Region Reflective and Nelder-Mead simplex are also available. For more details on the Levenberg-Marquardt algorithm, see the Wikipedia article.
This example function is provided as-is without any representation of accuracy.
Excel Usage
=LM_FIT(xdata, ydata, model_id, second_model, guess, bounds, lmfit_method)xdata(list[list], required): The xdata valueydata(list[list], required): The ydata valuemodel_id(str, required): The model_id valuesecond_model(str, optional, default: null): Optional second model to add for composite fitting (model name)guess(list[list], optional, default: null): Initial parameter guesses (single row with values for each parameter)bounds(list[list], optional, default: null): Parameter bounds (two rows - lower and upper bounds)lmfit_method(str, optional, default: “leastsq”): Fitting algorithm selection
Returns (list[list]): 2D list [param_names, fitted_values, std_errors], or error message string.
Examples
Example 1: Simple linear regression y = mx + b
Inputs:
| model_id | xdata | ydata |
|---|---|---|
| linear | 0 | 1 |
| 1 | 3 | |
| 2 | 5 | |
| 3 | 7 | |
| 4 | 9 |
Excel formula:
=LM_FIT("linear", {0;1;2;3;4}, {1;3;5;7;9})Expected output:
| slope | intercept |
|---|---|
| 2 | 1 |
| 0 | 0 |
Example 2: Constant model fit y = c
Inputs:
| model_id | xdata | ydata |
|---|---|---|
| constant | 0 | 5 |
| 1 | 5 | |
| 2 | 5 | |
| 3 | 5 |
Excel formula:
=LM_FIT("constant", {0;1;2;3}, {5;5;5;5})Expected output:
| c |
|---|
| 5 |
| 0 |
Example 3: Exponential decay model fit
Inputs:
| model_id | xdata | ydata |
|---|---|---|
| exponential | 0 | 10 |
| 1 | 3.679 | |
| 2 | 1.353 | |
| 3 | 0.498 | |
| 4 | 0.183 |
Excel formula:
=LM_FIT("exponential", {0;1;2;3;4}, {10;3.679;1.353;0.498;0.183})Expected output:
| amplitude | decay |
|---|---|
| 10 | 1 |
| 0.00026 | 0.000057 |
Example 4: Composite model with two Gaussians
Inputs:
| model_id | second_model | xdata | ydata | guess | |||||
|---|---|---|---|---|---|---|---|---|---|
| gaussian | gaussian | -3 | 0.054 | 0.4 | -1 | 1 | 0.4 | 1 | 1 |
| -2 | 0.246 | ||||||||
| -1 | 0.453 | ||||||||
| 0 | 0.484 | ||||||||
| 1 | 0.453 | ||||||||
| 2 | 0.246 | ||||||||
| 3 | 0.054 |
Excel formula:
=LM_FIT("gaussian", "gaussian", {-3;-2;-1;0;1;2;3}, {0.054;0.246;0.453;0.484;0.453;0.246;0.054}, {0.4,-1,1,0.4,1,1})Expected output:
| m1_amplitude | m1_center | m1_sigma | m2_amplitude | m2_center | m2_sigma |
|---|---|---|---|---|---|
| 1 | -1 | 1 | 1 | 1 | 1 |
| 0.0054 | 0.005 | 0.0024 | 0.0054 | 0.005 | 0.0024 |
Python Code
import micropip
await micropip.install(["lmfit"])
import numpy as np
import math
from lmfit import models as lmfit_models
def lm_fit(xdata, ydata, model_id, second_model=None, guess=None, bounds=None, lmfit_method='leastsq'):
"""
Fit data using lmfit's built-in models with optional model composition.
See: https://lmfit.github.io/lmfit-py/builtin_models.html
This example function is provided as-is without any representation of accuracy.
Args:
xdata (list[list]): The xdata value
ydata (list[list]): The ydata value
model_id (str): The model_id value Valid options: Gaussian, Lorentzian, Voigt, Pseudo-Voigt, Skewed Gaussian, Skewed Voigt, Pearson VII, Pearson IV, Moffat, Student's t, Exponential, Power Law, Quadratic, Polynomial 3rd Order, Polynomial 4th Order, Linear, Constant, Fermi, Step Function, Doniach, Sine Wave, Damped Oscillator, Damped Harmonic Oscillator, Lognormal, Exponential Gaussian, Breit-Wigner, Rectangle, Split Lorentzian.
second_model (str, optional): Optional second model to add for composite fitting (model name) Default is None.
guess (list[list], optional): Initial parameter guesses (single row with values for each parameter) Default is None.
bounds (list[list], optional): Parameter bounds (two rows - lower and upper bounds) Default is None.
lmfit_method (str, optional): Fitting algorithm selection Valid options: Levenberg-Marquardt, Trust Region Reflective, Nelder-Mead. Default is 'leastsq'.
Returns:
list[list]: 2D list [param_names, fitted_values, std_errors], or error message string.
"""
model_registry = {
'gaussian': {'class': lmfit_models.GaussianModel, 'param_hints': {'sigma': {'min': 0}}},
'lorentzian': {'class': lmfit_models.LorentzianModel, 'param_hints': {'sigma': {'min': 0}}},
'voigt': {'class': lmfit_models.VoigtModel, 'param_hints': {'sigma': {'min': 0}, 'gamma': {'min': 0}}},
'pseudovoigt': {'class': lmfit_models.PseudoVoigtModel, 'param_hints': {'sigma': {'min': 0}, 'fraction': {'min': 0, 'max': 1}}},
'skewedgaussian': {'class': lmfit_models.SkewedGaussianModel, 'param_hints': {'sigma': {'min': 0}}},
'skewedvoigt': {'class': lmfit_models.SkewedVoigtModel, 'param_hints': {'sigma': {'min': 0}, 'gamma': {'min': 0}}},
'pearson7': {'class': lmfit_models.Pearson7Model, 'param_hints': {'sigma': {'min': 0}, 'exponent': {'min': 0}}},
'pearson4': {'class': lmfit_models.Pearson4Model, 'param_hints': {'sigma': {'min': 0}, 'exponent': {'min': 0}}},
'moffat': {'class': lmfit_models.MoffatModel, 'param_hints': {'sigma': {'min': 0}, 'beta': {'min': 0}}},
'studentst': {'class': lmfit_models.StudentsTModel, 'param_hints': {'sigma': {'min': 0}}},
'exponential': {'class': lmfit_models.ExponentialModel, 'param_hints': {}},
'powerlaw': {'class': lmfit_models.PowerLawModel, 'param_hints': {}},
'quadratic': {'class': lmfit_models.QuadraticModel, 'param_hints': {}},
'polynomial3': {'class': lambda: lmfit_models.PolynomialModel(degree=3), 'param_hints': {}},
'polynomial4': {'class': lambda: lmfit_models.PolynomialModel(degree=4), 'param_hints': {}},
'linear': {'class': lmfit_models.LinearModel, 'param_hints': {}},
'constant': {'class': lmfit_models.ConstantModel, 'param_hints': {}},
'fermi': {'class': lmfit_models.FermiModel, 'param_hints': {'kt': {'min': 0}}},
'step': {'class': lmfit_models.StepModel, 'param_hints': {}},
'doniach': {'class': lmfit_models.DoniachModel, 'param_hints': {'sigma': {'min': 0}}},
'sine': {'class': lmfit_models.SineModel, 'param_hints': {}},
'dampedoscillator': {'class': lmfit_models.DampedOscillatorModel, 'param_hints': {}},
'dampedharmonicoscillator': {'class': lmfit_models.DampedHarmonicOscillatorModel, 'param_hints': {'sigma': {'min': 0}, 'gamma': {'min': 0}}},
'lognormal': {'class': lmfit_models.LognormalModel, 'param_hints': {'sigma': {'min': 0}}},
'exponentialgaussian': {'class': lmfit_models.ExponentialGaussianModel, 'param_hints': {'sigma': {'min': 0}}},
'breitwigner': {'class': lmfit_models.BreitWignerModel, 'param_hints': {'sigma': {'min': 0}, 'q': {'min': -10, 'max': 10}}},
'rectangle': {'class': lmfit_models.RectangleModel, 'param_hints': {'sigma1': {'min': 0}, 'sigma2': {'min': 0}}},
'splitlorentzian': {'class': lmfit_models.SplitLorentzianModel, 'param_hints': {'sigma': {'min': 0}, 'sigma_r': {'min': 0}}},
}
def normalize_to_2d_list(value):
"""Convert scalar [[x]] to 2D list for consistency with Excel."""
if not isinstance(value, list):
return [[value]]
return value
def validate_2d_list(data, name):
"""Validate and convert 2D list to 1D numpy array."""
if data is None:
return None
data = normalize_to_2d_list(data)
if not isinstance(data, list):
return None
if len(data) == 0:
return f"Error: {name} is empty."
first_col = [row[0] if isinstance(row, list) else row for row in data]
try:
arr = np.array(first_col, dtype=float)
# Check for NaN or inf in input data
if np.any(np.isnan(arr)) or np.any(np.isinf(arr)):
return f"Error: {name} contains NaN or infinite values."
return arr
except (ValueError, TypeError):
return f"Error: {name} contains non-numeric values."
try:
# Validate inputs
x_arr = validate_2d_list(xdata, "xdata")
y_arr = validate_2d_list(ydata, "ydata")
if isinstance(x_arr, str):
return x_arr
if isinstance(y_arr, str):
return y_arr
if len(x_arr) != len(y_arr):
return "Error: xdata and ydata must have same number of points."
if len(x_arr) < 3:
return "Error: Need at least 3 data points for fitting."
# Validate model IDs
if model_id not in model_registry:
valid = ', '.join(sorted(model_registry.keys()))
return f"Error: Invalid model_id '{model_id}'. Valid models: {valid}"
if second_model is not None and second_model not in model_registry:
valid = ', '.join(sorted(model_registry.keys()))
return f"Error: Invalid second_model '{second_model}'. Valid models: {valid}"
# Helper function to instantiate model from class or callable
def instantiate_model(model_class, prefix=''):
"""Instantiate a model with optional prefix."""
if callable(model_class):
if prefix:
# For lambdas like polynomial, we need to instantiate then get class
temp = model_class()
return temp.__class__(prefix=prefix)
return model_class()
return model_class(prefix=prefix)
# Get model info
model_info_1 = model_registry[model_id]
# Create composite model if second_model provided
if second_model is not None:
model_info_2 = model_registry[second_model]
# Create models with prefixes
model_1 = instantiate_model(model_info_1['class'], prefix='m1_')
model_2 = instantiate_model(model_info_2['class'], prefix='m2_')
model = model_1 + model_2
# Extract parameter names with prefixes
param_names = list(model_1.param_names) + list(model_2.param_names)
# Combine hints with prefixes
all_hints = {}
for p, hint in model_info_1['param_hints'].items():
all_hints['m1_' + p] = hint
for p, hint in model_info_2['param_hints'].items():
all_hints['m2_' + p] = hint
else:
# Single model without prefix
model = instantiate_model(model_info_1['class'])
param_names = list(model.param_names)
all_hints = model_info_1['param_hints']
# Initialize parameters
try:
params = model.guess(y_arr, x=x_arr)
except Exception:
# Fallback to make_params if guess fails
params = model.make_params()
# Apply user-provided guesses
if guess is not None:
guess_normalized = normalize_to_2d_list(guess)
if not isinstance(guess_normalized, list) or len(guess_normalized) == 0:
return "Error: guess must be a non-empty list."
if len(guess_normalized[0]) != len(param_names):
return f"Error: guess has {len(guess_normalized[0])} values, expected {len(param_names)}."
for i, pname in enumerate(param_names):
try:
val = float(guess_normalized[0][i])
if math.isnan(val) or math.isinf(val):
return f"Error: guess contains invalid value (NaN or inf) for parameter '{pname}'."
params[pname].value = val
except (KeyError, ValueError, TypeError):
return f"Error: Could not set guess for parameter '{pname}'."
# Apply bounds
if bounds is not None:
bounds_normalized = normalize_to_2d_list(bounds)
if not isinstance(bounds_normalized, list) or len(bounds_normalized) != 2:
return "Error: bounds must have 2 rows (lower and upper)."
lower_bounds = bounds_normalized[0]
upper_bounds = bounds_normalized[1]
if not isinstance(lower_bounds, list) or not isinstance(upper_bounds, list):
return "Error: bounds rows must be lists."
if len(lower_bounds) != len(param_names) or len(upper_bounds) != len(param_names):
return f"Error: bounds have wrong dimensions. Expected {len(param_names)} columns."
for i, pname in enumerate(param_names):
try:
lb = float(lower_bounds[i]) if lower_bounds[i] is not None else -np.inf
ub = float(upper_bounds[i]) if upper_bounds[i] is not None else np.inf
if not math.isnan(lb) and not math.isnan(ub) and lb > ub:
return f"Error: Lower bound exceeds upper bound for parameter '{pname}'."
params[pname].min = lb
params[pname].max = ub
except (KeyError, ValueError, TypeError):
return f"Error: Could not set bounds for parameter '{pname}'."
# Apply parameter hints (positivity constraints, etc.)
for pname, hints in all_hints.items():
if pname in params:
for key, value in hints.items():
setattr(params[pname], key, value)
# Perform fit
if lmfit_method is None:
result = model.fit(y_arr, params, x=x_arr)
else:
result = model.fit(y_arr, params, x=x_arr, method=lmfit_method)
# Extract results
fitted_vals = []
std_errors = []
for pname in param_names:
if pname not in result.params:
return f"Error: Parameter '{pname}' missing from fit results."
val = float(result.params[pname].value)
stderr = result.params[pname].stderr
stderr = float(stderr) if stderr is not None else np.nan
# Validate no NaN/Inf in fitted values
if math.isnan(val) or math.isinf(val):
return "Error: Fitting produced invalid numeric values (NaN or inf)."
fitted_vals.append(val)
std_errors.append(stderr)
return [param_names, fitted_vals, std_errors]
except Exception as e:
return f"Error: {str(e)}"