Fitting Gaussian curve on top of polynomial curve?

Hi I want to fit a Gaussian curve on top of a polynomail curve. I am playing around with the Curve Fitting toolbox and one of the tutorials they give is this image:

As I understand you can put in a custom equation but I'm not sure what sort of equation to use . If the larger decaying curve (the one whose peak is beyond 0 on the X-axis) is the one I want to be polynomial, and they the center hump at ~125 is the Gaussian curve, how do I go about doing this? It's OK if I need to do multiple steps in the code.

Answer 1

[Saarthak Gupta]

Hi Tiffany,

It seems like you want to fit a combination of a Gaussian curve and a polynomial curve to your data.

The problem you are referring to is known as segmented regression or piecewise regression, which involves finding a piecewise continuous function to best describe the data.

Typically, to solve a segmented regression problem, the least squares method is applied separately to each segment.

The two regression lines are made to fit the data set as closely as possible while minimizing the sum of squares of the differences (SSD) between observed (y) and calculated (Yr) values of the dependent variable.

Since the exact data is not provided, I have prepared a model dataset by combining samples from a polynomial and a Gaussian curve:

Refer to the following code:

% Polynomial samples (straight line)
x1 = randi(62, [1 100]);    % random integer samples in (0,62]
y1 = (6/-119000).*(x1-120);
et1 = -1+2.*rand(size(y1)); % error term/noise
y1 = y1+et1./7000;

% Gaussian samples
x2 = randi(250,[1 500]);    % random integer samples in (0,250]
x2(x2<61) = [];
y2 = normpdf(x2,125,50);    % calculate density of samples
et2 = -1+2.*rand(size(y2)); % error term/noise
y2 = y2+et2./7000;

% Combine samples
x = [x1 x2];
y = [y1 y2];

scatter(x,y,".k")

fitSegPoly(x1,y1);  % Fit polynomial segment
fitSegGauss(x2,y2); % Fit Gaussian segment

fitSegPoly.m

function [fitresult, gof] = fitSegPoly(x1, y1)
%CREATEFIT(X1,Y1)
%  Create a fit.
%
%  Data for 'polySegmentFit' fit:
%      X Input: x1
%      Y Output: y1
%  Output:
%      fitresult : a fit object representing the fit.
%      gof : structure with goodness-of fit info.
%
%  See also FIT, CFIT, SFIT.

%  Auto-generated by MATLAB on 19-Dec-2023 13:48:39


%% Fit: 'polySegmentFit'.
[xData, yData] = prepareCurveData( x1, y1 );

% Set up fittype and options.
ft = fittype( 'poly1' );

% Fit model to data.
[fitresult, gof] = fit( xData, yData, ft );

% Plot fit with data.
figure( 'Name', 'polySegmentFit' );
h = plot( fitresult, xData, yData );
legend( h, 'y1 vs. x1', 'polySegmentFit', 'Location', 'NorthEast', 'Interpreter', 'none' );
% Label axes
xlabel( 'x1', 'Interpreter', 'none' );
ylabel( 'y1', 'Interpreter', 'none' );
grid on

fitSegGauss.m

function [fitresult, gof] = fitSegGauss(x2, y2)
%CREATEFIT(X2,Y2)
%  Create a fit.
%
%  Data for 'gaussSegmentFit' fit:
%      X Input: x2
%      Y Output: y2
%  Output:
%      fitresult : a fit object representing the fit.
%      gof : structure with goodness-of fit info.
%
%  See also FIT, CFIT, SFIT.

%  Auto-generated by MATLAB on 19-Dec-2023 13:49:54


%% Fit: 'gaussSegmentFit'.
[xData, yData] = prepareCurveData( x2, y2 );

% Set up fittype and options.
ft = fittype( 'gauss1' );
opts = fitoptions( 'Method', 'NonlinearLeastSquares' );
opts.Display = 'Off';
opts.Lower = [-Inf -Inf 0];
opts.StartPoint = [0.00810733661388935 125 32.2342882092758];

% Fit model to data.
[fitresult, gof] = fit( xData, yData, ft, opts );

% Plot fit with data.
figure( 'Name', 'gaussSegmentFit' );
h = plot( fitresult, xData, yData );
legend( h, 'y2 vs. x2', 'gaussSegmentFit', 'Location', 'NorthEast', 'Interpreter', 'none' );
% Label axes
xlabel( 'x2', 'Interpreter', 'none' );
ylabel( 'y2', 'Interpreter', 'none' );
grid on

The following curves are produced for the segments:

Once you identify a suitable breakpoint, you can fit the respective segments using a polynomial and a Gaussian curve separately. Alternatively, you can also try to fit a higher order Gaussian curve to your data, without segmenting.

Refer to the following MATLAB documentation for further reference:

• https://www.mathworks.com/help/curvefit/fit.html