Fit x and y data

10 views (last 30 days)
Thor
Thor on 22 Feb 2021
Edited: Thor on 14 May 2021
Hello everyone,
I got this data and I want to create a script that plots a gaussian curve that fits the major peak only. I have seen many examples.However, they involve superimposing normal distribution funcitons or are for multiple peaks , so the fitting i get isnt correct.

Accepted Answer

Star Strider
Star Strider on 22 Feb 2021
Try this:
D1 = readmatrix('Test1.xls');
x = D1(:,1);
y = D1(:,2);
gausfcn = @(b,x) b(1).*exp(-(x-b(2)).^2/b(3));
[maxy,idx] = max(y);
Lv = y >= 1.5; % Restrict Region Of Fit To Region Of Symmetry In ‘y’
B = fminsearch(@(b)norm(y(Lv) - gausfcn(b,x(Lv))), [maxy, x(idx), 1/12.5]);
figure
plot(x, y)
hold on
plot(x, gausfcn(B,x), '-r')
hold off
grid
text(x(idx), 0.4, sprintf('$y = %5.2f \\cdot e^{\\frac{-(x-%7.2f)^2}{%7.2f}}$', B), 'Interpreter','latex', 'HorizontalAlignment','center', 'FontSize',12)
producing:
.
  10 Comments
Thor
Thor on 16 Mar 2021
Hi @Star Strider, I was referring to this
'' If you have the Global Optimization Toolbox, I will post the code I used for this so you can use the ga function to experiment with it ''
Star Strider
Star Strider on 16 Mar 2021
Sure!
I ran it again to be certain that it works.
D1 = readmatrix('Test1.xls');
x = D1(:,1);
y = D1(:,2);
[maxy,idx] = max(y);
Lv = y >= 1.08; % Restrict Region Of Fit To Region Of Symmetry In ‘y’
gausfcn2a = @(b,x) b(1).*exp(-(x-b(2)).^2/b(3)) + b(4).*exp(-(x-b(5)).^2/b(6));
% gausfcn2m = @(b,x) b(1).*exp(-(x-b(2)).^2/b(3)) .* b(4).*exp(-(x-b(5)).^2/b(6));
ftns = @(b)norm(y(Lv) - gausfcn2a(b,x(Lv)))
PopSz = 50;
Parms = 6;
opts = optimoptions('ga', 'PopulationSize',PopSz, 'InitialPopulationMatrix',randi(1E+4,PopSz,Parms)*1E-3+[maxy, x(idx), 1/12.5, maxy, x(idx), rand], 'MaxGenerations',2E5, 'PlotFcn',@gaplotbestf, 'PlotInterval',1);
t0 = clock;
fprintf('\nStart Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t0)
[B,fval,exitflag,output] = ga(ftns, Parms, [],[],[],[],zeros(Parms,1),Inf(Parms,1),[],[],opts)
t1 = clock;
fprintf('\nStop Time: %4d-%02d-%02d %02d:%02d:%07.4f\n', t1)
GA_Time = etime(t1,t0)
DT_GA_Time = datetime([0 0 0 0 0 GA_Time], 'Format','HH:mm:ss.SSS');
fprintf('\nElapsed Time: %23.15E\t\t%s\n\n', GA_Time, DT_GA_Time)
fprintf(1,'\tRate Constants:\n')
for k1 = 1:length(B)
fprintf(1, '\t\tB(%d) = %8.5f\n', k1, B(k1))
end
figure
plot(x, y)
hold on
plot(x, gausfcn2a(B,x), '-r')
hold off
grid
% text(x(idx), 0.4, sprintf('$y = %5.2f \\cdot e^{\\frac{-(x-%7.2f)^2}{%7.2f}}$', B), 'Interpreter','latex', 'HorizontalAlignment','center', 'FontSize',12)
% text(x(idx), 0.2, sprintf('$y = %5.2f \\cdot e^{(\\frac{-(x-%7.2f)}{%7.2f})^2}$', B(1:2),sqrt(B(3))), 'Interpreter','latex', 'HorizontalAlignment','center', 'FontSize',12)
Experiment with it to get the result you want. It may be necessary to tweak the tolerances to produce a slightly better result than this code produces. You can use any of the objective function variations in the ‘ftns’ function, or create your own, to get a suitable result.
Adding a bit of documentation —
PopSz = population size (rows of the 'InitialPopulationMatrix')
Parms = Number of Parameters in the objective function (columns of the 'InitialPopulationMatrix')
This vector: ‘[maxy, x(idx), 1/12.5, maxy, x(idx), rand]’ scales the columns of the 'InitialPopulationMatrix' to different magnitudes. That makes it easier for the ga algorithm to converge. Change that (or eliminate it entirely) as necessary.
If you have other questions, post them and I will do my best to provide appropriate replies.

Sign in to comment.

More Answers (0)

Community Treasure Hunt

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

Start Hunting!