Fitting Gaussian to a curve with multiple peaks

696 views (last 30 days)
I have a curve with a few peaks (the figure is attached here). I want to fit a few Gaussian to it. How can I do that?
  1 Comment
dpb on 21 Jul 2018
There are pre-built models in the Curve Fitting Toolbox if you have it...

Sign in to comment.

Answers (4)

Image Analyst
Image Analyst on 12 Oct 2020
Attached is a demo for how to fit any specified number of Gaussians to noisy data. Here is an example where I created a signal from 6 component Gaussians by summing then, and then added noise to the summed curve. The input data is the dashed line (upper most curve), and the Gaussians it thought would sum to fit it best are shown in the solid color curves below the dashed curve.
Here is the main part of the program:
% Run optimization
[parameter, fval, flag, output] = fminsearch(@(lambda)(fitgauss(lambda, tFit, y)), startingGuesses, options);
The rest is just setup/initialization and plotting code.
Thor on 22 Feb 2021
Edited: Thor on 22 Feb 2021
Hi @Image Analyst, I got a similar issue. I got some UV data and i want to fit a gaussian to the major peak. I have attempted to use your code, but the fitting is not as nice as the one your have shown above. Please , could you provide me with some guidance? I would be happy to open a question. Many thanks.

Sign in to comment.

Image Analyst
Image Analyst on 21 Jul 2018
Gaussian Mixture Models
Gaussian mixture models (GMM) are composed of k multivariate normal density components, where k is a positive integer. Each component has a d-dimensional mean (d is a positive integer), d-by-d covariance matrix, and a mixing proportion. Mixing proportion j determines the proportion of the population composed by component j, j = 1,...,k.
You can fit a GMM using the Statistics and Machine Learning Toolbox™ function fitgmdist by specifying k and by supplying X, an n-by-d matrix of data. The columns of X correspond to predictors, features, or attributes, and the rows correspond to observations or examples. By default, fitgmdist fits full covariance matrices that are different among components (or unshared).
Rik on 12 Oct 2020
Comment posted as flag by qi lu:
Hi, this demo seems perfect. Could you please tell me how to fit such equation. I use 'fit' in the curve fitting box, the fitting result is horrible.

Sign in to comment.

Image Analyst
Image Analyst on 24 Jul 2018
Edited: Image Analyst on 25 Nov 2020
One way to do it is to set up the sum of two Gaussians with an offset and a linear ramp. Then you can use fitnlm, with your best guesses as to the parameters. See code:
% Uses fitnlm() to fit a non-linear model (sum of two gaussians on a ramp) through noisy data.
% Requires the Statistics and Machine Learning Toolbox, which is where fitnlm() is contained.
% Initialization steps.
clc; % Clear the command window.
close all; % Close all figures (except those of imtool.)
clear; % Erase all existing variables. Or clearvars if you want.
workspace; % Make sure the workspace panel is showing.
format long g;
format compact;
fontSize = 20;
% Create the X coordinates from 0 to 20 every 0.5 units.
X = linspace(0, 80, 400);
mu1 = 10; % Mean, center of Gaussian.
sigma1 = 3; % Standard deviation.
mu2 = 50; % Mean, center of Gaussian.
sigma2 = 9; % Standard deviation.
% Define function that the X values obey.
a = 6 % Arbitrary sample values I picked.
b = 35
c = 25
m = 0.1
Y = a + m * X + b * exp(-(X - mu1) .^ 2 / sigma1)+ c * exp(-(X - mu2) .^ 2 / sigma2); % Get a vector. No noise in this Y yet.
% Add noise to Y.
Y = Y + 1.2 * randn(1, length(Y));
% Now we have noisy training data that we can send to fitnlm().
% Plot the noisy initial data.
plot(X, Y, 'b.', 'LineWidth', 2, 'MarkerSize', 15);
grid on;
% Convert X and Y into a table, which is the form fitnlm() likes the input data to be in.
tbl = table(X', Y');
% Define the model as Y = a + exp(-b*x)
% Note how this "x" of modelfun is related to big X and big Y.
% x((:, 1) is actually X and x(:, 2) is actually Y - the first and second columns of the table.
modelfun = @(b,x) b(1) + b(2) * x + b(3) * exp(-(x(:, 1) - b(4)).^2/b(5)) + b(6) * exp(-(x(:, 1) - b(7)).^2/b(8));
beta0 = [6, 0.1, 35, 10, 3, 25, 50, 9]; % Guess values to start with. Just make your best guess.
% Now the next line is where the actual model computation is done.
mdl = fitnlm(tbl, modelfun, beta0);
% Now the model creation is done and the coefficients have been determined.
% YAY!!!!
% Extract the coefficient values from the the model object.
% The actual coefficients are in the "Estimate" column of the "Coefficients" table that's part of the mode.
coefficients = mdl.Coefficients{:, 'Estimate'}
% Let's do a fit, but let's get more points on the fit, beyond just the widely spaced training points,
% so that we'll get a much smoother curve.
X = linspace(min(X), max(X), 1920); % Let's use 1920 points, which will fit across an HDTV screen about one sample per pixel.
% Create smoothed/regressed data using the model:
yFitted = coefficients(1) + coefficients(2) * X + coefficients(3) * exp(-(X - coefficients(4)).^2 / coefficients(5)) + coefficients(6) * exp(-(X - coefficients(7)).^2 / coefficients(8));
% Now we're done and we can plot the smooth model as a red line going through the noisy blue markers.
hold on;
plot(X, yFitted, 'r-', 'LineWidth', 2);
grid on;
title('Exponential Regression with fitnlm()', 'FontSize', fontSize);
xlabel('X', 'FontSize', fontSize);
ylabel('Y', 'FontSize', fontSize);
legendHandle = legend('Noisy Y', 'Fitted Y', 'Location', 'northeast');
legendHandle.FontSize = 25;
% Set up figure properties:
% Enlarge figure to full screen.
set(gcf, 'Units', 'Normalized', 'OuterPosition', [0 0 1 1]);
% Get rid of tool bar and pulldown menus that are along top of figure.
% set(gcf, 'Toolbar', 'none', 'Menu', 'none');
% Give a name to the title bar.
set(gcf, 'Name', 'Demo by ImageAnalyst', 'NumberTitle', 'Off')
Image Analyst
Image Analyst on 27 Feb 2021
I haven't done it in 2-D. You'd have to experiment around with it and try to adapt it to 2-D. You can do that as well as I can.

Sign in to comment.

Star Strider
Star Strider on 25 Jul 2018
Another option: Area under each peak (link).
There seems to be a ‘hidden’ peak at about x=50, so you might have to isolate it by first subtracting the others, or include it manually, since I doubt findpeaks will detect it as a peak.
Also, you may need to determine if the baseline variations are real (a true wandering baseline). or are the result of some peaks having wide distributions. Since this could be an iterative process, it may be best to ignore the baseline variations until you are convinced a wandering baseline exists.

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!