Bi-exponential curve fitting using fminsearch

Hello all,

I am having issues while trying to modify a monoexponential curve fitting routine to perform a biexponential fit. I am fairly new to matlab so this is probably the root of my issue. The curve fitting seems to fit the data well, but the output parameters seem incorrect. The exponential fit formula is below:

FittedCurve =A*(exp(-lambda*xdata)) + B*(exp(-lambda2*xdata));

where A and B are the amplitudes each exponential, and lambda and lambda2 are the rate constants for each. These parameters are changing when fitting the same data multiple times...

Here is the curve fitting function code I am using:

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

function [estimates, model] = curvefit_biexp(xdata, ydata)

start_point = rand(1,4);

model =@efun; options = optimset('Display','off','TolFun',1e-16,'TolX',1e-16);

estimates = fminsearch(model, start_point,options);

% expfun accepts curve parameters as inputs, and outputs sse,

% and the FittedCurve.

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

end

Also here is a set of sample data:

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

Any help with this issue would be greatly appreciated!!! Thanks in advance...

[Terence]

Another little detail: This curve fitting function is going to be used for analyzing metabolic data from phosphorus magnetic resonance spectroscopy for those familiar. The idea to make measurements of metabolic capacity. Here is a similar equation used previously in research:

PCr and pH recoveries were fitted as a sum of exponentials with amplitudes A1 and A2, and rate constants k1 and k2. PCr kinetics were fitted as:

y = [Rest − Depl] [A1(1 − e^−tk1) +A2(1 − e^−tk2)] + Depl,

where Rest is the fully recovered PCr level and Depl is the initial PCr level.

This is what I am hoping to mimick with the above function...

Thanks All!

Terence

Answer 1

[Matt J]

Your model is rather ill-conditioned. For example, because the role of (A,lambda) is interchangeable with (B,lambda2), you have non-unique solutions. As another example, if your ydata happens to be a single exponential

ydata=C*exp(-L*xdata)

then you will have a continuum of non-unique solutions consisting of lambda=lambda2=L and any A,B satisfying A+B=C.

You need to add physical info to the model so that the roles of the 2 terms are more distinguishable in some way.

Also, your objective function has large flat regions where fminsearch can get stuck if your lambdas get too large. Since your xdata are mostly much greater than 10, for example, your true lambdas must be much less than 1/10 Otherwise, your fitted curve will be essentially zero almost everywhere. Initializing with rand() will not put you in that zone reliably. You might even need to consider doing a constrained optimization to keep the lambdas from getting unrealistically large.