Azzera filtri
Azzera filtri

How to fit multiple gaussian in a curve ?

149 visualizzazioni (ultimi 30 giorni)
Hello everyone,
I want to get individual gaussian functions to get the area under each peak
So, how I can fit a multiple gaussian in a curve with multiple peaks

Risposta accettata

John D'Errico
John D'Errico il 25 Feb 2023
Modificato: John D'Errico il 25 Feb 2023
Easy enough. You NEED good starting values though. Crappy starting values --> crappy results.
load dsp.mat
plot(f,pxx)
grid on
So there is one mode around x==100, a second mode near x==450, a third mode near x==600, a 4th mode near 750. What matters most is to have good starting values for the modes.
mdl = fittype('gauss4')
mdl =
General model Gauss4: mdl(a1,b1,c1,a2,b2,c2,...,a4,b4,c4,x) = a1*exp(-((x-b1)/c1)^2) + a2*exp(-((x-b2)/c2)^2) + a3*exp(-((x-b3)/c3)^2) + a4*exp(-((x-b4)/c4)^2)
fittedmdl = fit(f,pxx,mdl,'start',[1e-5 110 20 1e-5 450 100 1e-5 600 50 1e-5 750 50])
fittedmdl =
General model Gauss4: fittedmdl(x) = a1*exp(-((x-b1)/c1)^2) + a2*exp(-((x-b2)/c2)^2) + a3*exp(-((x-b3)/c3)^2) + a4*exp(-((x-b4)/c4)^2) Coefficients (with 95% confidence bounds): a1 = 0.0006381 (0.0005812, 0.0006951) b1 = 110 (108.5, 111.5) c1 = 20 (17.94, 22.06) a2 = 0.0002224 (0.0001959, 0.0002489) b2 = 450 (435.6, 464.4) c2 = 100 (78.58, 121.4) a3 = 9.539e-05 (5.252e-05, 0.0001383) b3 = 600 (579.4, 620.6) c3 = 50 (19.8, 80.2) a4 = 5.798e-05 (2.164e-05, 9.432e-05) b4 = 750 (724, 776) c4 = 50 (11.86, 88.14)
plot(fittedmdl,f,pxx)
And, as you can see, despite my using decent starting values so it found 4 gaussian modes where I told it, they all fit poorly. This points out the next issue. YOUR CURVE IS NOT FIT WELL BY A SUM OF GAUSSIAN MODES! They look vaguely gaussian, but they are not in fact accurately fit by gaussians. I suppose I could have added another mode out around 900, and maybe split that mode near x==450 into two modes. Even then, look at the first pek. It is pretty much completely alone, yet a Gaussian does not fit well. Even so, I am sure you will want to try harder.
mdl = fittype('gauss6')
mdl =
General model Gauss6: mdl(a1,b1,c1,a2,b2,c2,...,a6,b6,c6,x) = a1*exp(-((x-b1)/c1)^2) + a2*exp(-((x-b2)/c2)^2) + a3*exp(-((x-b3)/c3)^2) + a4*exp(-((x-b4)/c4)^2) + a5*exp(-((x-b5)/c5)^2) + a6*exp(-((x-b6)/c6)^2)
fittedmdl = fit(f,pxx,mdl,'start',[1e-5 110 20 1e-5 400 100 1e-5 450 100 1e-5 600 50 1e-5 750 50 1e-5 900 50])
fittedmdl =
General model Gauss6: fittedmdl(x) = a1*exp(-((x-b1)/c1)^2) + a2*exp(-((x-b2)/c2)^2) + a3*exp(-((x-b3)/c3)^2) + a4*exp(-((x-b4)/c4)^2) + a5*exp(-((x-b5)/c5)^2) + a6*exp(-((x-b6)/c6)^2) Coefficients (with 95% confidence bounds): a1 = 0.0006381 (0.0005833, 0.0006929) b1 = 110 (108.6, 111.4) c1 = 20 (18.01, 21.99) a2 = 8.221e-05 (-0.128, 0.1282) b2 = 400 (-3.778e+04, 3.858e+04) c2 = 100 (-5884, 6084) a3 = 0.0001487 (-0.126, 0.1263) b3 = 450 (-2.298e+04, 2.388e+04) c3 = 100 (-4664, 4864) a4 = 0.0001066 (-0.0001708, 0.000384) b4 = 600 (563, 637) c4 = 50 (-9.948, 109.9) a5 = 5.774e-05 (2.224e-05, 9.324e-05) b5 = 750 (723.9, 776.1) c5 = 50 (10.42, 89.58) a6 = 1.613e-05 (-1.886e-05, 5.111e-05) b6 = 900 (810.3, 989.7) c6 = 50 (-81.77, 181.8)
plot(fittedmdl,f,pxx)
Still, total, complete, unambiguous crap. Note that between the first and second peak, there is a wide tail, far too wide for a Gaussian to fit. And the first peak itself has wider tails than would be suggested by a Gaussian term.
Using a sum of Gaussians is a poor solution to solve this problem. I'm sorry, but it is. (I'm not going to suggest a good solution, because that itself is not remotely trivial.)
Just because something has a bump does not mean it is well represented by a Gaussian curve shape. The problem is, you will not get a good estimate of the area under those peaks. And where there are several things happening at once, things are far more difficult to decide what area corresponds to each vaguely gaussian "bump".
  4 Commenti
Abdelhakim Souidi
Abdelhakim Souidi il 25 Feb 2023
thank you @John D'Errico, if I want to only plot for example the first gaussian
a1*exp(-((x-b1)/c1)^2)
How can I do that
John D'Errico
John D'Errico il 26 Feb 2023
Was that not my example in my last comment? That is, I plotted the first Gaussian. I did this:
Gfun = @(x,a,b,c) a*exp(-((x-b)/c).^2);
Do you understand this is the Gaussian shaped function you are using? Next, what did I do?
fplot(@(x) Gfun(x,fmdl.a1,fmdl.b1,fmdl.c1),[0,200])
I plotted the FIRST gaussian. Why are you asking how to do exactly what I just showed you how to do?

Accedi per commentare.

Più risposte (2)

Image Analyst
Image Analyst il 26 Feb 2023
It looks like your data has about 6 peaks. See my attached demo, that just happens to use 6 peaks but that is a variable number you can change in the code to whatever you want. But you should really know what the theory says. Like are you expecting two peaks or 6? Use whatever the theory says you should be expecting. I have attached another demo that uses fitnlm to fit exactly two Gaussians.
In the above plot you can see the upper 2 curves match pretty well. One is the original data and the other (dashed one) is the sum of all 6 Gaussians.
  1 Commento
Abdelhakim Souidi
Abdelhakim Souidi il 26 Feb 2023
Thank you @Image Analyst
I already have your code but still incapable to adapt it to my data, more precisely I don't know where to modify
could you please specify how to adapt it to my data

Accedi per commentare.


Alex Sha
Alex Sha il 16 Apr 2023
For the summation of 6 Gaussians function:
Sum Squared Error (SSE): 2.61218021364296E-9
Root of Mean Square Error (RMSE): 4.49993989076388E-6
Correlation Coef. (R): 0.999145442992317
R-Square: 0.998291616252313
Parameter Best Estimate
--------- -------------
a1 0.000616781092742577
a2 0.000218333421684105
a3 -6.3828311170045E-6
a4 -0.000138418678090533
a5 0.000138852036452085
a6 0.000133702730236499
b1 105.976485712389
b2 106.829210687263
b3 7.3291143621259E-14
b4 309.94152870094
b5 379.758692434269
b6 425.217299664439
c1 -9.29831066559539
c2 25.1779734434923
c3 520.573845470695
c4 115.976664646099
c5 -119.557220394127
c6 295.98790157476

Prodotti


Release

R2022a

Community Treasure Hunt

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

Start Hunting!

Translated by