How to do Curve Fitting with a Piecewise Gaussian Model

Hi,

I am trying to fit a piecewise Gaussian function to an intensity profile (blue curve in the picture) and I'm running into Problems extracting the best fit parameters. This is my first time trying to do a Curve Fit so I decided to ask you for some advice.

I implemented the piecewise Model:

    function [ value ] = PiecewiseGaussianModel( A1,A2,m1,m2,sigma1,sigma2,I1,I2,A,B,x )
    value=zeros(size(x));
    for i=1:length(x)
        if x(i) < A || x(i) > B
            value(i)=-A1*exp(-((x(i)-m1)^2)/(2*(sigma1)^2))+I1;
        elseif A <= x(i) && x(i) <= B
            value(i)= A2*exp(-((x(i)-m2)^2)/(2*(sigma2)^2))+I2;
        end
    end
    end

and the minimization problem:

    function [ error ] = GaussianError( parameter,x )
    load('profile2.mat'); %y intensity Data
    x=1:1:numel(y);
    x=x';
    [ A,B ] = FindBoundaries( y );
    fx=PiecewiseGaussianModel(parameter(1),parameter(2),parameter(3),parameter(4), ...
                              parameter(5),parameter(6),parameter(7),parameter(8),A,B,x); 
    fx=fx';
    error=fx-y;
    end;

    function [ A,B ] = FindBoundaries( y )
    [~, ind1]=min(y);
    copy_y=y;
    copy_y(ind1)=inf;
    [~, ind2]=min(copy_y);
    if ind1 > ind2
       A=ind2;
       B=ind1;
    else
       A=ind1;
       B=ind2;
    end
    end

profile2.mat returns y and this is how y' looks like:

y'=

Columns 1 through 13

   130   120   106    95    90    84    85    90    94   101   112   122   121

Columns 14 through 26

   106    96    90    87    82    81    83    88    97   110   124   136   147

Columns 27 through 28

   153   151

I tried solving this problem by using the Matlab functions lsqcurvefit and lsqnonlin, but I don't get any good results. Is this problem solvable with these functions? And if so, is there a good solution getting the start points for the fitting?

Thanks in advance

Answer 1

[John D'Errico]

If there were good ways to find good starting values for some general nonlinear optimization, then the code would do it for you.

Is the problem solvable with those functions? I think you may have a serious problem, in that the breakpoints create a discontinuous function. It is not clear if you are trying to use an optimizer to move the breakpoints. So as a breakpoint is moved by the code, points can move from one Gaussian segment to another, across that discontinuity. So your sum of squares objective is in turn non-smooth, and in fact, discontinuous itself.

By the way, it looks like you are passing in the parameters as a list of separate arguments. Note that lsqcurvefit and lsqnonlin require one vector as an argument. That vector would pack all of your arguments together. Of course, I don't know how you called the codes, so how can I know what you actually did?

Should you be fitting this curve as a set of pieces of Gaussians? Shiver. I don't think you should do so. I don't know why you are trying to use Gaussian segments, except that you got the idea that you could do so. Your boss or a colleague suggested the idea, or you read it somewhere. These things never work well. You never get very good information from that little amount of data, and a Gaussian never has quite the right curve shape. Here you have maybe a half dozen or so points that fall into each segment, so not enough to fit a Gaussian very well anyway.

A better idea might be to fit a spline. For example, using my SLM toolbox on the data you posted, this seems reasonable:

slm = slmengine(0:27,y,'plot','on','knots',15);

Now that we can see your data, I would point out that the individual peaks do not have terribly Gaussian-like shapes. So any fits that expect them to be so will be fraught with problems.

Find my SLM toolbox on the File Exchange.

If your end goal was simply a code that will give you a nice curve fit, then you can stop there.

If your end goal was to extract information about the shape of the peaks, like "standard deviations" of those peaks as a means of defining their widths, then you still need to do some work, but at least now you have a continuous function to work with.