CHROMA_PEAKS
Overview
The CHROMA_PEAKS function fits specialized peak shape models to chromatographic data using nonlinear least-squares optimization. Chromatography is an analytical technique used to separate mixtures, and the resulting output typically shows peaks corresponding to different compounds. Accurate characterization of peak shape is essential for quantitative analysis, peak deconvolution, and compound identification.
This function implements seven distinct peak models commonly used in chromatographic analysis:
- Beta Distribution Asymmetric Peak: A flexible model using beta distribution parameters to capture asymmetric peak shapes with controllable skewness.
- Bigaussian Asymmetric Peak: A piecewise Gaussian model with different widths (
w1,w2) on each side of the peak center, ideal for simple asymmetric peaks. - Chesler-Cram Exponential Tail Peak: Combines a Gaussian core with an exponential tailing function to model peaks with significant trailing edges.
- Edgeworth-Cramer Skewed Peak: Extends the Gaussian distribution using Edgeworth series expansion to incorporate skewness and kurtosis corrections.
- Exponentially Modified Gaussian (EMG) Peak: A convolution of a Gaussian and an exponential distribution, widely used as a theoretical model for chromatographic peak shapes. See the Wikipedia article on EMG distributions for mathematical details.
- Giddings Bessel Dispersion Peak: Based on Bessel functions to model dispersion effects in chromatographic systems, using the modified Bessel function of the first kind.
- Weibull Amplitude Asymmetric Peak: Uses the Weibull distribution to model asymmetric peaks with flexible tail behavior.
The underlying optimization uses scipy.optimize.curve_fit from the SciPy library, which implements the Levenberg-Marquardt algorithm for unconstrained problems or trust-region reflective algorithms when parameter bounds are specified. Each model includes intelligent initial parameter guesses derived from the input data and appropriate bounds to ensure physically meaningful results. The function returns fitted parameter values along with standard errors estimated from the covariance matrix.
This example function is provided as-is without any representation of accuracy.
Excel Usage
=CHROMA_PEAKS(xdata, ydata, chroma_peaks_model)xdata(list[list], required): The xdata valueydata(list[list], required): The ydata valuechroma_peaks_model(str, required): The chroma_peaks_model value
Returns (list[list]): 2D list [param_names, fitted_values, std_errors], or error string.
Examples
Example 1: Demo case 1
Inputs:
| chroma_peaks_model | xdata | ydata |
|---|---|---|
| beta_distribution_asymmetric_peak | 1.2 | 0.86160756 |
| 1.6 | 1.19043718 | |
| 2 | 1.3 | |
| 2.4 | 1.20340707 | |
| 2.8 | 0.95652922 |
Excel formula:
=CHROMA_PEAKS("beta_distribution_asymmetric_peak", {1.2;1.6;2;2.4;2.8}, {0.86160756;1.19043718;1.3;1.20340707;0.95652922})Expected output:
"non-error"
Example 2: Demo case 2
Inputs:
| chroma_peaks_model | xdata | ydata |
|---|---|---|
| bigaussian_asymmetric_peak | 0.01 | 1.3331896701095287 |
| 2.0075 | 5.013805785170162 | |
| 4.005 | 0.5789147596326245 | |
| 6.0024999999999995 | 0.01 | |
| 8 | 0.03795354352729467 |
Excel formula:
=CHROMA_PEAKS("bigaussian_asymmetric_peak", {0.01;2.0075;4.005;6.0024999999999995;8}, {1.3331896701095287;5.013805785170162;0.5789147596326245;0.01;0.03795354352729467})Expected output:
"non-error"
Example 3: Demo case 3
Inputs:
| chroma_peaks_model | xdata | ydata |
|---|---|---|
| chesler_cram_exponential_tail_peak | 0.01 | 0.6881970247706027 |
| 2.0075 | 3.476308512165623 | |
| 4.005 | 0.2972799078474512 | |
| 6.0024999999999995 | 0.01 | |
| 8 | 0.026603606259349277 |
Excel formula:
=CHROMA_PEAKS("chesler_cram_exponential_tail_peak", {0.01;2.0075;4.005;6.0024999999999995;8}, {0.6881970247706027;3.476308512165623;0.2972799078474512;0.01;0.026603606259349277})Expected output:
"non-error"
Example 4: Demo case 4
Inputs:
| chroma_peaks_model | xdata | ydata |
|---|---|---|
| edgeworth_cramer_skewed_peak | 0.01 | 0.593946685207516 |
| 2.0075 | -0.3106033211311715 | |
| 4.005 | -0.06975777024010682 | |
| 6.0024999999999995 | 0.03426042007874261 | |
| 8 | 0.006244911955196433 |
Excel formula:
=CHROMA_PEAKS("edgeworth_cramer_skewed_peak", {0.01;2.0075;4.005;6.0024999999999995;8}, {0.593946685207516;-0.3106033211311715;-0.06975777024010682;0.03426042007874261;0.006244911955196433})Expected output:
"non-error"
Example 5: Demo case 5
Inputs:
| chroma_peaks_model | xdata | ydata |
|---|---|---|
| exponentially_modified_gaussian_peak | 0.01 | 1 |
| 2.0075 | 1 | |
| 4.005 | 1 | |
| 6.0024999999999995 | 1 | |
| 8 | 1 |
Excel formula:
=CHROMA_PEAKS("exponentially_modified_gaussian_peak", {0.01;2.0075;4.005;6.0024999999999995;8}, {1;1;1;1;1})Expected output:
"non-error"
Example 6: Demo case 6
Inputs:
| chroma_peaks_model | xdata | ydata |
|---|---|---|
| giddings_bessel_dispersion_peak | 0.01 | 1 |
| 2.0075 | 1 | |
| 4.005 | 1 | |
| 6.0024999999999995 | 1 | |
| 8 | 1 |
Excel formula:
=CHROMA_PEAKS("giddings_bessel_dispersion_peak", {0.01;2.0075;4.005;6.0024999999999995;8}, {1;1;1;1;1})Expected output:
"non-error"
Example 7: Demo case 7
Inputs:
| chroma_peaks_model | xdata | ydata |
|---|---|---|
| weibull_amplitude_asymmetric_peak | 0.01 | 0.08678315093623744 |
| 2.0075 | 1.006439786471626 | |
| 4.005 | 0.3007064729763167 | |
| 6.0024999999999995 | 0.03250067264273272 | |
| 8 | 0.014649108934044 |
Excel formula:
=CHROMA_PEAKS("weibull_amplitude_asymmetric_peak", {0.01;2.0075;4.005;6.0024999999999995;8}, {0.08678315093623744;1.006439786471626;0.3007064729763167;0.03250067264273272;0.014649108934044})Expected output:
"non-error"
Python Code
import numpy as np
from scipy import special as sc
from scipy.optimize import curve_fit as scipy_curve_fit
import math
def chroma_peaks(xdata, ydata, chroma_peaks_model):
"""
Fits chroma_peaks models to data using scipy.optimize.curve_fit. See https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.html for details.
See: https://docs.scipy.org/doc/scipy/reference/generated/scipy.optimize.curve_fit.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
chroma_peaks_model (str): The chroma_peaks_model value Valid options: Beta Distribution Asymmetric Peak, Bigaussian Asymmetric Peak, Chesler Cram Exponential Tail Peak, Edgeworth Cramer Skewed Peak, Exponentially Modified Gaussian Peak, Giddings Bessel Dispersion Peak, Weibull Amplitude Asymmetric Peak.
Returns:
list[list]: 2D list [param_names, fitted_values, std_errors], or error string.
"""
def _validate_data(xdata, ydata):
"""Validate and convert both xdata and ydata to numpy arrays."""
for name, arg in [("xdata", xdata), ("ydata", ydata)]:
if not isinstance(arg, list) or len(arg) < 2:
raise ValueError(f"{name}: must be a 2D list with at least two rows")
vals = []
for i, row in enumerate(arg):
if not isinstance(row, list) or len(row) == 0:
raise ValueError(f"{name} row {i}: must be a non-empty list")
try:
vals.append(float(row[0]))
except Exception:
raise ValueError(f"{name} row {i}: non-numeric value")
if name == "xdata":
x_arr = np.asarray(vals, dtype=np.float64)
else:
y_arr = np.asarray(vals, dtype=np.float64)
if x_arr.shape[0] != y_arr.shape[0]:
raise ValueError("xdata and ydata must have the same number of rows")
return x_arr, y_arr
# Model definitions dictionary
models = {
'beta_distribution_asymmetric_peak': {
'params': ['y0', 'xc', 'A', 'w1', 'w2', 'w3'],
'model': lambda x, y0, xc, A, w1, w2, w3: y0 + A * np.power(1.0 + ((w2 + w3 - 2.0) / (w2 - 1.0)) * ((x - xc) / w1), w2 - 1.0) * np.power(1.0 - ((w2 + w3 - 2.0) / (w3 - 1.0)) * ((x - xc) / w1), w3 - 1.0),
'guess': lambda xa, ya: (float(np.min(ya)), float(np.median(xa)), float(np.max(ya)), float(max(np.max(xa) - np.min(xa), 1e-3)), 2.0, 2.0),
'bounds': ([-np.inf, -np.inf, 0.0, 0.0, 1.0000001, 1.0000001], np.inf),
},
'bigaussian_asymmetric_peak': {
'params': ['y0', 'xc', 'H', 'w1', 'w2'],
'model': lambda x, y0, xc, H, w1, w2: y0 + H * np.exp(-0.5 * np.square((x - xc) / np.where(x < xc, w1, w2))),
'guess': lambda xa, ya: (float(np.min(ya)), float(np.median(xa)), float(np.ptp(ya) if np.ptp(ya) else 1.0), float(max((np.max(xa) - np.min(xa)) / 6.0, 1e-3)), float(max((np.max(xa) - np.min(xa)) / 6.0, 1e-3))),
'bounds': ([-np.inf, -np.inf, -np.inf, 0.0, 0.0], np.inf),
},
'chesler_cram_exponential_tail_peak': {
'params': ['y0', 'xc1', 'A', 'w', 'k2', 'xc2', 'B', 'k3', 'xc3'],
'model': lambda x, y0, xc1, A, w, k2, xc2, B, k3, xc3: y0 + A * (np.exp(-np.square(x - xc1) / (2.0 * w)) + B * (1.0 - 0.5 * (1.0 - np.tanh(k2 * (x - xc2)))) * np.exp(-0.5 * k3 * (np.abs(x - xc3) + (x - xc3)))),
'guess': lambda xa, ya: (float(np.min(ya)), float(np.median(xa)), float(np.ptp(ya) if np.ptp(ya) else 1.0), float(max(np.square(np.max(xa) - np.min(xa)) / 8.0, 1e-3)), 1.0, float(np.median(xa)), 0.5, 1.0, float(np.median(xa))),
'bounds': ([-np.inf, -np.inf, -np.inf, 0.0, 0.0, -np.inf, -np.inf, 0.0, -np.inf], np.inf),
},
'edgeworth_cramer_skewed_peak': {
'params': ['y0', 'xc', 'A', 'w', 'a3', 'a4'],
'model': lambda x, y0, xc, A, w, a3, a4: y0 + (A / (w * np.sqrt(2.0 * np.pi))) * np.exp(-0.5 * np.square((x - xc) / w)) * (1.0 + (a3 / 6.0) * (x - xc) / w * ((np.square((x - xc) / w)) - 3.0) + (a4 / 24.0) * (np.power((x - xc) / w, 4) - 6.0 * np.square((x - xc) / w) + 3.0) + (10.0 * np.square(a3) / 720.0) * (np.power((x - xc) / w, 6) - 15.0 * np.power((x - xc) / w, 4) + 45.0 * np.square((x - xc) / w) - 15.0)),
'guess': lambda xa, ya: (float(np.min(ya)), float(np.median(xa)), float(np.ptp(ya) if np.ptp(ya) else 1.0), float(max((np.max(xa) - np.min(xa)) / 6.0, 1e-3)), 0.0, 0.0),
'bounds': ([-np.inf, -np.inf, 0.0, 0.0, -np.inf, -np.inf], np.inf),
},
'exponentially_modified_gaussian_peak': {
'params': ['y0', 'A', 'xc', 'w', 't0'],
'model': lambda x, y0, A, xc, w, t0: y0 + (A / t0) * np.exp(0.5 * np.square(w / t0) - (x - xc) / t0) * 0.5 * (1.0 + sc.erf(((x - xc) / w - w / (2.0 * t0)) / np.sqrt(0.5))),
'guess': lambda xa, ya: (float(np.min(ya)), float(np.ptp(ya) if np.ptp(ya) else 1.0), float(np.median(xa)), float(max((np.max(xa) - np.min(xa)) / 6.0, 1e-3)), 1.0),
'bounds': ([-np.inf, 0.0, -np.inf, 0.0, 0.0], np.inf),
},
'giddings_bessel_dispersion_peak': {
'params': ['y0', 'xc', 'w', 'A'],
'model': lambda x, y0, xc, w, A: y0 + A / w * np.sqrt(np.clip(xc, 1e-9, None) / np.clip(x, 1e-9, None)) * sc.iv(1, 2.0 * np.sqrt(np.clip(xc, 1e-9, None) * np.clip(x, 1e-9, None)) / w) * np.exp(-(x + xc) / w),
'guess': lambda xa, ya: (float(np.min(ya)), float(np.median(xa)), float(max((np.max(xa) - np.min(xa)) / 4.0, 1e-3)), float(np.ptp(ya) if np.ptp(ya) else 1.0)),
'bounds': ([-np.inf, -np.inf, 0.0, 0.0], np.inf),
},
'weibull_amplitude_asymmetric_peak': {
'params': ['y0', 'xc', 'A', 'w1', 'w2'],
'model': lambda x, y0, xc, A, w1, w2: y0 + A * np.exp(-np.exp(-(x - xc) / w2) - (x - xc) / w1),
'guess': lambda xa, ya: (float(np.min(ya)), float(np.median(xa)), float(np.ptp(ya) if np.ptp(ya) else 1.0), float(max((np.max(xa) - np.min(xa)) / 4.0, 1e-3)), 1.0),
'bounds': ([-np.inf, -np.inf, -np.inf, 0.0, 0.0], np.inf),
}
}
# Validate model parameter
if chroma_peaks_model not in models:
return f"Invalid model: {str(chroma_peaks_model)}. Valid models are: {', '.join(models.keys())}"
model_info = models[chroma_peaks_model]
# Validate and convert input data
try:
x_arr, y_arr = _validate_data(xdata, ydata)
except ValueError as e:
return f"Invalid input: {e}"
# Perform curve fitting
try:
p0 = model_info['guess'](x_arr, y_arr)
bounds = model_info.get('bounds', (-np.inf, np.inf))
if bounds == (-np.inf, np.inf):
popt, pcov = scipy_curve_fit(model_info['model'], x_arr, y_arr, p0=p0, maxfev=10000)
else:
popt, pcov = scipy_curve_fit(model_info['model'], x_arr, y_arr, p0=p0, bounds=bounds, maxfev=10000)
fitted_vals = [float(v) for v in popt]
for v in fitted_vals:
if math.isnan(v) or math.isinf(v):
return "Fitting produced invalid numeric values (NaN or inf)."
except ValueError as e:
return f"Initial guess error: {e}"
except Exception as e:
return f"curve_fit error: {e}"
# Calculate standard errors
std_errors = None
try:
if pcov is not None and np.isfinite(pcov).all():
std_errors = [float(v) for v in np.sqrt(np.diag(pcov))]
except Exception:
pass
return [model_info['params'], fitted_vals, std_errors] if std_errors else [model_info['params'], fitted_vals]