How to Fit Gaussian Peaks: Mathematics and Solver Initialization

Scientific data parsing often requires resolving peak profiles from raw sensor readings. Whether analyzing thermal Doppler broadening, optical spectra, or chromatography signals, the Gaussian profile is the most common mathematical model used to describe peak shapes.

The Mathematics of a Gaussian Peak

A single Gaussian peak in experimental datasets is defined by the following equation:

Where:

• A is the peak amplitude (height above baseline).
• μ (mu) is the center position of the peak.
• σ (sigma) is the standard deviation (controlling width).
• y₀ is a constant vertical offset (baseline offset).

For chemists and physicists, the peak width is usually expressed as FWHM (Full Width at Half Maximum). The mathematical relationship between FWHM and σ is:

Why Solver Initialization is Critical

To fit this equation to noisy experimental data, we use non-linear least squares solvers (such as the Levenberg-Marquardt algorithm). However, non-linear optimization requires initial guesses for the parameters. If the initial guesses are poor, the solver might diverge, fail to fit, or get trapped in local minima.

Here is a robust heuristic algorithm to automatically initialize the parameters:

1. Baseline (y₀): Take the median or average of the 5% lowest data points.
2. Center (μ): Locate the x-coordinate of the maximum y-value in the dataset.
3. Amplitude (A): Maximum y-value minus the estimated baseline y₀.
4. Width (σ): Estimate the FWHM by finding the distance between points where the intensity drops to half of the peak amplitude, then divide by 2.3548.

Implementing it in Python (SciPy)

If you are using AltaiPlot's built-in Python Scripting Console, or working in a traditional Jupyter Notebook, here is how you can implement this robust heuristic:

import numpy as np
from scipy.optimize import curve_fit

# 1. Define the Gaussian model
def gaussian(x, A, mu, sigma, y0):
    return A * np.exp(-((x - mu) ** 2) / (2 * sigma ** 2)) + y0

x_data = np.linspace(-10, 10, 200)
y_noisy = gaussian(x_data, A=5.0, mu=1.2, sigma=1.5, y0=1.0) + np.random.normal(0, 0.2, len(x_data))

# 2. Heuristic Initialization (Crucial for convergence!)
guess_y0 = np.mean(np.sort(y_noisy)[:int(len(y_noisy)*0.05)]) # Lowest 5% mean
max_idx = np.argmax(y_noisy)
guess_mu = x_data[max_idx]
guess_A = y_noisy[max_idx] - guess_y0

# Estimate FWHM
half_max = guess_y0 + guess_A / 2
above_half = np.where(y_noisy > half_max)[0]
fwhm_guess = (x_data[above_half[-1]] - x_data[above_half[0]]) if len(above_half) > 1 else 2.0
guess_sigma = fwhm_guess / 2.35482

initial_guesses = [guess_A, guess_mu, guess_sigma, guess_y0]

# 3. Fit the model
popt, pcov = curve_fit(gaussian, x_data, y_noisy, p0=initial_guesses)
print(f"Fitted Parameters:\nA: {popt[0]:.3f}, mu: {popt[1]:.3f}, sigma: {popt[2]:.3f}, y0: {popt[3]:.3f}")

This script robustly extracts bounds, calculates means for the baseline, and finds the argmax for the center. If you execute this in AltaiPlot's Python console, you can easily append the result popt directly back to your active plot!

The Point-and-Click Alternative (Zero-Code)

While scripting is incredibly powerful for custom batch jobs, writing boilerplate initialization code for every single routine analysis is a bottleneck. For those who prefer a faster workflow, AltaiPlot provides an interactive Fit Panel that automates this entire script:

• Point-and-Click Anchoring: Simply drag a selection box over the peak region. AltaiPlot's engine automatically executes the parameter guess heuristic in the background using prominence-based peak detection.
• CPU-based Levenberg-Marquardt Solver: The non-linear solver runs on your local machine with bound constraints and numeric Jacobian caching, fitting multiple peaks simultaneously.
• Visual Error & Residuals: Residual curves are plotted in real-time below the chart, letting you see exactly how well the Gaussian model matches your raw data.
• Speed Modes: Choose between Quick (5K subsampled), Fit (100K cap), or Advanced (streaming LM for 100M+ points) depending on your dataset size.