How to fit multiple gaussian in a curve ?

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

Answer 1

[John D'Errico]

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".

[Abdelhakim Souidi]

Thank you for your answer, I appreciate it. Could you just tell me how to know the subgaussiansthat participate to form this curve

[John D'Errico]

Whatever floats your boat. You can extract the coefficients one at a time, like this:

fittedmdl.a1

or extract the entire set of coefficients as

coeffvalues(fittedmdl)

All I can do is tell you the area estimates you will get will be wholly inaccurate. And the problem will not be because fit does not produce a sufficiently accurate fit. The problem will more likely be because you will be using the wrong model for that fit, and at best, that your starting values need to be amazingly good.

You would get a far better estimate for the area under that first peak, for example, just by using trapz.

load dsp.mat
plot(f,cumtrapz(f,pxx))
grid on
xline(200,'r')

So the first peak lives below x==200. cumtrapz tells me the area in that peak is roughly 0.025.

Again, this is a difficult thing to do, because you are trying to solve what falls into a class of ill-posed problems.

Ok. if you REALLY think the fit is a good thing to do...

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)

fmdl = fit(f,pxx,mdl,'start',[1e-5 110 20 1e-5 450 100 1e-5 600 50 1e-5 750 50])

fmdl = 
     General model Gauss4:
     fmdl(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)

Gfun = @(x,a,b,c) a*exp(-((x-b)/c).^2);
fplot(@(x) Gfun(x,fmdl.a1,fmdl.b1,fmdl.c1),[0,200])

Surely you can compute the area under that mode. The location of the mode is irrelevant. So the area under a peak is simply:

syms x a c real
int(a*exp(-(x/c)^2),[-inf,inf])

ans = 

[Abdelhakim Souidi]

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]

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?

Answer 2

[Image Analyst]

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.

[Abdelhakim Souidi]

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

Answer 3

[Alex Sha]

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 

[James]

Can you link your code please?

[Walter Roberson]

@James

@Alex Sha used the third party toolbox named "1stOpt" http://www.7d-soft.com/en/ . Alex has a hobby of running requested fits through that software and posting the results. I have always been amazed at how well that toolbox fits data, and how quickly it operates, I have considered buying a license myself... unfortunately I just cannot justify the price.

[Sam Chak]

Hi @James

You should zoom in on the fitted curve to evaluate the quality of the fit from a human perspective. Even with a very high R-squared value, there can be deviations between the curve and the data points.

More importantly, you should investigate why the built-in fit() function did not produce a similar result when the 1stOpt results were assigned as initial fitting values to the coefficients of the gauss6 model.

As @John D'Errico stated, 'Crappy starting values, crappy results.' Does this not conversely imply that 'Quality in, quality out' (QIQO)?

format long g
load('dsp.mat');

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      ;

%% Fit Alex's gauss6 model to the data
x   = linspace(f(1), f(end), 100001);
y   = 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);

figure
plot(f, pxx, '.'), hold on
plot(x, y), grid on, hold off
legend('Data', 'Fitted Model', 'FontSize', 14)
title('Fitting result by 1stOpt')
xlabel('f'), ylabel('pxx')

figure
tl  = tiledlayout(2, 2);
nexttile
plot(f, pxx, '.'), hold on
plot(x, y), grid on, hold off
xlim([0 250])

nexttile
plot(f, pxx, '.'), hold on
plot(x, y), grid on, hold off
xlim([250 500])

nexttile
plot(f, pxx, '.'), hold on
plot(x, y), grid on, hold off
xlim([500 750])

nexttile
plot(f, pxx, '.'), hold on
plot(x, y), grid on, hold off
xlim([750 1000])

xlabel(tl, 'f'), ylabel(tl, 'pxx')

%% Fit six Gaussian functions using Alex's best values as starting points
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)

[Gauss6, gof] = fit(f, pxx, mdl, 'StartPoint', [a1 b1 c1 a2 b2 c2 a3 b3 c3 a4 b4 c4 a5 b5 c5 a6 b6 c6])

Gauss6 = 
     General model Gauss6:
     Gauss6(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
       b1 =         106
       c1 =         0.5
       a2 =   0.0005616
       b2 =       106.8
       c2 =       25.18
       a3 =  -6.225e-05
       b3 =   7.329e-14
       c3 =       520.6
       a4 =  -4.072e-05
       b4 =       309.9
       c4 =         116
       a5 =           0
       b5 =       379.8
       c5 =         0.5
       a6 =   0.0001972  (-0.002226, 0.002621)
       b6 =       425.2  (-1031, 1881)
       c6 =         296  (-891.9, 1484)

gof = struct with fields:
           sse: 2.50125236632583e-07
       rsquare: 0.836403305203325
           dfe: 111
    adjrsquare: 0.81134795554978
          rmse: 4.74697852344357e-05

figure
plot(Gauss6, f, pxx), grid on
title('MATLAB Fit using 1stOpt results as initial values')