Exponential Curve fitting

Hello,

I am fairly new to Matlab and have been teaching myself for a few months. I have written a code to curve fit some data and calculate time and rate constants for the exponential recovery for some data. I have been running into some problems curve fitting the data, and I cannot figure out where the problem is. The code will fit some data, and have issues with others. Any help would be great! Below, I will post the curve fitting function code, and also the data I am trying to fit.

Thanks in Advance!!!!

Here is the Curve Fitting Function:

function [estimates, model] = curvefit(xdata, ydata)
% fits data to the curve y(x)=A-B*e(-lambda*x))
start_point = [1,1,1];

 model =@efun;
 estimates = fminsearch(model, start_point);

% expfun accepts curve parameters as inputs, and outputs sse,
% the sum of squares error for A -B* exp(-lambda * xdata) - ydata,
% and the FittedCurve. 

      function [sse,FittedCurve] = efun(v)
          A=v(1);
          B=v(2);
          lambda=v(3);
          FittedCurve =A - B*exp(-lambda*xdata);
          ErrorVector=FittedCurve-ydata;
          sse = sum(ErrorVector .^2);
      end
  end

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

Here is some data I am trying to fit with this code

    xdata=[0;10.7;21.2;31.7;41.9;52.2;62.6;72.6;82.9;92.9;103.2;113.2;123.5;133.8;144;154.3;164.9;175.2;185.5;196;206.9;217.2;227.4;237.7;248.3];

    ydata=[-2.5808;-2.1815;-1.6986;-1.2984;-0.9516;-0.7136;-0.4505;-0.3751;-0.3165;-0.2603;-0.2263;-0.1996;-0.1658;-0.1983;-0.2140;-0.2365;-0.1680;-0.1751;-0.1757;-0.1827;-0.1968;-0.1811;-0.1550;-0.1901;-0.1510];

[the cyclist]

Can you be more specific about what you mean by "have issues with other [data]"? MATLAB throws an error? The fit is not as good as you think it should be? Be as specific as you can.

Answer 1

[Matt Tearle]

Looks a whole lot like a typical local-vs-global minimum issue. The minimum that fminsearch finds depends heavily on the initial guess. Change the line start_point = [1,1,1]; to start_point = rand(1,3); then run the following lines a few times:

xdata=[0;10.7;21.2;31.7;41.9;52.2;62.6;72.6;82.9;92.9;103.2;113.2;123.5;133.8;144;154.3;164.9;175.2;185.5;196;206.9;217.2;227.4;237.7;248.3];
ydata=[-2.5808;-2.1815;-1.6986;-1.2984;-0.9516;-0.7136;-0.4505;-0.3751;-0.3165;-0.2603;-0.2263;-0.1996;-0.1658;-0.1983;-0.2140;-0.2365;-0.1680;-0.1751;-0.1757;-0.1827;-0.1968;-0.1811;-0.1550;-0.1901;-0.1510];

[c,f] = curvefit(xdata,ydata)
[~,ymod] = f(c);
plot(xdata,ydata,'o',xdata,ymod)

Sometimes you'll see the weird result, sometimes it will work nicely.

So... what can you do?

1. find some way to get an estimate of the coefficients (for start_point)
2. use a global optimization routine (if you have Global Opt Toolbox), such as multistart
3. fake multistart. Use a random starting point, then do something like:

xdata=[0;10.7;21.2;31.7;41.9;52.2;62.6;72.6;82.9;92.9;103.2;113.2;123.5;133.8;144;154.3;164.9;175.2;185.5;196;206.9;217.2;227.4;237.7;248.3];
ydata=[-2.5808;-2.1815;-1.6986;-1.2984;-0.9516;-0.7136;-0.4505;-0.3751;-0.3165;-0.2603;-0.2263;-0.1996;-0.1658;-0.1983;-0.2140;-0.2365;-0.1680;-0.1751;-0.1757;-0.1827;-0.1968;-0.1811;-0.1550;-0.1901;-0.1510];
numattempts = 50;

err = Inf;
for k=1:numattempts
  [c,f] = curvefit(xdata,ydata);
  [thiserr,thismodel] = f(c);
  if thiserr<err
      err = thiserr;
      coeffs = c;
      ymodel = thismodel;
  end
end

disp(coeffs)
plot(xdata,ydata,'o',xdata,ymodel)

[Terence Ryan]

Hello Matt,

Thanks for your help. I think I understand what you are saying, but let me just repeat it for my sake. The curvefit I have written searches for local minima, which can be altered heavily by the starting_point. I have thought about provide initial 'guesses' for these coefficients, but I thought this method could bias my results in way that is not ideal...

The code you have provided above uses a random starting_point. The added code the create a 'fake multistart' seems to work great.

Thanks for the Help! I think this has solved my problems!!! If anyone has other suggestions send them along.

[Matt Tearle]

Exactly! You can also incorporate things like The Cyclist suggested. Note that start_point = rand(1,3); makes initial values in the range [0,1], so you might want to scale these, to allow a greater range.

Answer 2

[the cyclist]

I think a more stringent tolerance for the fit will get you what you want. Use these lines in place of your call to fminsearch:

 options = optimset('Display','iter','TolFun',1e-16);
 estimates = fminsearch(model, start_point,options);

[Terence Ryan]

Sorry I should have been more clear about the 'issues' I am having. For some data the fit does not appear to be accurate. Sometimes, the data look exponential, but the curve fit code returns a square curve (Time constant = 0.29 sec). Other issues are just with the accuracy of the curve. The residuals of the curve fit are all substantially positive numbers towards the end of the data (i.e.- the curve never touches the baseline data.) Your addition for more stringent tolerance looks like it has helped a great deal, so thanks for the help! Any more advice would be so helpful.

I am also getting the following output after running the curve fit:

Exiting: Maximum number of function evaluations has been exceeded

- increase MaxFunEvals option.

Current function value: 0.069213

Any suggestions on this?