MATLAB: Exponential Curve fitting

Question

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];

Answer 1

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)