How to fit a Gauss curve in the tails of a curve

I need to fit a Gauss curve to the tails of curve as shown in the picture, to calculate the area below the curve. I tried to use the code below without success. I'm using the wrong method ?

There is a better way to do it ?

I successfully did it for a Gaussian curve using the examples in the page https://www.mathworks.com/help/curvefit/fitoptions.html, shown in the picture. However, when I applied the same method to my curve it didn't work.

I tried using both, the std from the original curve and the std from a fitted Gaussian curve.

Here is my code (I attached the gdata.m file):

load('gdata.mat');

x = (1:1:129)';

% mystd = std(gdata)

gfit1= fit(x,gdata,'gauss2');

YHat = gfit1(x);

mystd = std(YHat)

gfit2= fit(x,gdata,'gauss2','Exclude', gdata > (2*mystd));

[dpb]

Not sure which area you're speaking of -- I just fit a single gaussian to the left/right tails from the intersection points where the inflection appears discontinuous that I identified by trial and error as L=39 and R=89 and got the following baseline approximation

but not sure what you're looking for next.

Perhaps some background on what generated the raw data and the end objective would help...

[Sergio Soares]

Perfect, this is what I want. However, I have 150 curves to find the Gaussian fit and the trial and error method for each curve will take a long time. I'm looking for an automatic way of doing it. This is why I tried the two standard deviations approach (after fit a Gaussian curve to this curve).

I need to calculate the area under the gaussian curve to find the scatter value of a PET measurement (phantom). This curve is the sum of the lines of a 2D sinogram.

Answer 1

[dpb]

If the data typically have such a discontinuity at the tails, a technique we used in the old online "SulfurMeter" which was a gamma-spec device for measuring elemental S in coal for power plant emissions control was something similar to the following illustration --

First a replotting of your data for a look-see more closely...

Note the sharpness of the discontinuity is enhanced by the first derivative. Occasionally, the second will be even more helpful but here it seems the inevitable additional noise is more than the enhancement and the location appears the same anyway...so, now what? Let's try fitting a breakpoint to the lefthand section--

>> [mxL,estL]=max(d);       % left tail max location, magnitude
>> cest=[0 0 estL-5 mxL]    % estimate for fitting breakpoint by piecewise linear fit
cest =
         0           0          37        8546
>> coeffL=nlinfit(1:estL,d(1:estL),@piecewise,cest)  % do the fit
coeffL =
 1.0e+03 *
 -0.0914    0.0123    0.0381    2.2369
>> [coeffL(3) round(coeffL(3))]   % the third coefficient is the breakpoint location
ans =
 38.0916   38.0000
>>

How well did it work?

>> figure
>> plot(d(1:Lest))
>> hold on
>> plot(1:Lest,piecewise(coeffL,1:Lest),'r-')
>> legend('1st Diff','Piecewise','location','northwest')

Looks pretty good. You can do the same thing from the left except use min and work from the RH end to the front in selecting the RH locations. I haven't done it here but you'll probably want to "add one" to the returned location of the breakpoint for use with the raw data since the difference vector is one element shorter--or, you can augment/prepend the difference vector with the mean of the first few elements or similar before fitting.

NB: The initial estimates are

0,0    -- first section intercept, slope
estL-5 -- breakpoint at arbitrary point a few positions off the max
mxL    -- slope for second section

The breakpoint location estimate may need some fine-tuning on setting the distance from the located maximum; it does need to be at least "reasonably close" to the knee you're trying to locate to get a good fit.

The slope is an approximation to the slope of the second line; it would be approximately (max-mean)/(xmax-startx) if were to try to actually compute; I just assumed only a width of 1 point and zero mean so it was (max-0)/(1). As you'll note, the fitted line slope is roughly 1/4 that so dividing by the offset number may make solution a little quicker but typically the solution is not terribly sensitive to this initial guess if the location is reasonable.

Another "trick" we used in the S-meter was to do an iterative fitting of subtracting a baseline then refitting the residuals.

piecewise is my utility function for the purpose--

function y=piecewise(coef,x)
% return piecewise linear fit given a1,b1,c,b2 as coefficient array and vector x
% c is breakpoint, a1 is intercept of first section, b1,b2 are two segment slopes

a1=coef(1); b1=coef(2);  % reduce parentheses clutter....
c=coef(3); b2=coef(4);
ix=x>c;                  % location > breakpoint
y(ix)=[a1+c*(b1-b2)+b2*x(ix)];
y(~ix)=[a1+b1*x(~ix)];
y=y(:);
>>

[Sergio Soares]

Hello, after some coffee I did the right breakpoint too and found 89. This are very good results without using try and error. I'm sorry for the lame mistakes and questions I did.

Thanks a lot for the help. Your solution to find the breakpoints is easy to implement and produce excellent results.

This is my final code:

load('sinogram_curves_27cm.mat');
gdata=squeeze(out(:,1,1));
temp=diff(gdata);
[mxL,estL]=max(temp);       % left tail max location, magnitude
cest=[0 0 estL-5 mxL];    % estimate for fitting breakpoint by piecewise linear fit
coeffL=nlinfit(1:estL,temp(1:estL),@piecewise,cest);  % do the fit
[coeffL(3) round(coeffL(3))]   % the third coefficient is the breakpoint location
[mxR,estR]=min(temp);
cest1=[0 128 estR+5 mxR];
coeffR=nlinfit(estR:128,temp(estR:end),@piecewise,cest1);  % do the fit
[coeffR(3) round(coeffR(3))]   % the third coefficient is the breakpoint location

[dpb]

Will be interesting to see how well it holds up for a larger population of data. As noted, the basic idea is embedded in the code for the online sulfurmeter and was a significant improvement in identifying the onset of "real" peak from the background. The SM was particularly painful in that it was gamma-spec w/ NaI crystal so there was a lot of overlap from the inherent width in peak resolution. Unfortunately, the cost for the required cooling at the time precluded Ge in the industrial setting.

Yours seems much better behaved from at least this spectrum...

PS. Sorry for the typo in changing variables in midstream...I just pasted from the command window and didn't check some later code was consistent withe earlier had written before...

I cleaned up the questions caused mostly from that and another false start or two since you'd already found the solutions on your own and didn't really contribute anything in the end...

Answer 2

[Image Analyst]

Just get a new gdata where you exclude the parts you don't want and include the parts you want. Like

x = [gx(1:100), gx(300:end)];
y = [gdata(1:100), gdata(300:end)];

Attach your x,y data if you need help with that.

[Sergio Soares]

Thanks for the help.

My data was attached. The file gdata.mat has the y values and the x values are in the code on the question.

x = (1:1:129)';

I did as you told me and worked well. However, I have 150 curves to calculate. It is a lot of work to see each one and extract the parts. This is why I was trying to use 2 standard deviations in the gfit function.

[Image Analyst]

You can find the middle part if you define some kind of rule as to what is the middle part. Like it's above 5000 or something. Or you might be able to use findchangepts(), or use triangle thresholding (attached).

[Sergio Soares]

Thanks for the help again.

I tried the findchangepts(), but it didn't work well at first. I had to adjust by trial and error the value of the maxnumofchanges to 10. Even this way it only found 41 and 89 for the break points which gave me a wrong fit. However, using the same thinking of my original code I got a good fit in the end, by excluding the values higher than 2 STD (peak).

In the pictures you can see the first try (with a peak) and the final result of the gauss fitting curve.

I think this can work. I will try it to the other curves and see what happens.

This is my final code for that:

load('sinogram_curves_27cm.mat');

gdata=squeeze(out(:,1,1));

xaxis=findchangepts(gdata,'MaxNumChanges',10)

temp = gdata([1:xaxis(1),xaxis(10):end]);

mystd = std(temp)

x = (1:1:size(temp,1))';

gfit = fit(x,temp,'gauss2');

gfit1= fit(x,temp,'gauss2','Exclude', temp > (2*mystd));

YHat = gfit(x);

YHat1 = gfit1(x);

figure; subplot(2,1,1) plot (x,YHat);

subplot(2,1,2) plot (x,YHat1);

[Sergio Soares]

I forgot to show the scatter plot of the first try.

[dpb]

Did you not try the empirical fitting to determine the breakpoint location I illustrated below? If not, why not?

[Sergio Soares]

I'm trying that now.