MATLAB - Exponential Curve Fitting without Toolbox

Question

I have data points (x, y) that I need to fit an exponential function to,

y = A + B * exp(C * x),

but I can neither use the Curve Fitting Toolbox nor the Optimization Toolbox.

User rayryeng was good enough to help me with working code:

x = [0    0.0036    0.0071    0.0107    0.0143    0.0178    0.0214    0.0250    0.0285    0.0321    0.0357    0.0392    0.0428    0.0464    0.0464];
y = [1.3985    1.3310    1.2741    1.2175    1.1694    1.1213    1.0804    1.0395    1.0043    0.9691    0.9385    0.9080    0.8809    0.7856    0.7856];

M = [ones(numel(x),1), x(:)]; %// Ensure x is a column vector
lny = log(y(:)); %// Ensure y is a column vector and take ln

X = M\lny; %// Solve for parameters
A = exp(X(1)); %// Solve for A
b = X(2); %// Get b

xval = linspace(min(x), max(x));
yval = A*exp(b*xval);
plot(x,y,'r.',xval,yval,'b');

However, this code only fits the equation without offset

y = A * exp(B * x).

How can I extend this code to fit the three-parameter equation?

---

In another attempt, I managed to fit the function using fminsearch:

function [xval, yval] = curve_fitting_exponential_1_optimized(x,y,xval)

start_point = rand(1, 3);
model = @expfun;
est = fminsearch(model, start_point);

    function [sse, FittedCurve] = expfun(params)
        A = params(1);
        B = params(2);
        C = params(3);
        FittedCurve = A + B .* exp(-C * x);
        ErrorVector = FittedCurve - y;
        sse = sum(ErrorVector .^ 2);
    end

yval = est(1)+est(2) * exp(-est(3) * xval);
end

The problem here is that the result depends on the starting point which is randomly chosen, so I don't get a stable solution. But since I need the function for automatization, I need something stable. How can I get a stable solution?

Answer 1

How to adapt rayryeng's code for three parameters?

rayryeng used the strategy to linearize a nonlinear equation so that standard regression methods can be applied. See also Jubobs' answer to a similar question.

This strategy does no longer work if there is a non-zero offset A. We can fix the situation by getting a rough estimate of the offset. As rubenvb mentioned in the comments, we could estimate A by min(y), but then the logarithm gets applied to a zero. Instead, we could leave a bit of space between our guess of A and the minimum of the data, say half its range. Then we subtract A from the data and use rayreng's method:

x = x(:);  % bring the data into the standard, more
y = y(:);  % convenient format of column vectors
Aguess = min(y) - (max(y) - min(y) / 2;
guess = [ones(size(x)), -x] \ log(y - Aguess);
Bguess = exp(guess(1));
Cguess = guess(2);

For the given data, this results in

Aguess =  0.4792
Bguess =  0.9440
Cguess = 21.7609

Other than for the two-parameter situation, we cannot expect this to be a good fit. Its SSE is 0.007331.

How to get a stable solution?

This guess is however useful as a starting point for the nonlinear optimization:

start_point = [Aguess, Bguess, Cguess];
est = fminsearch(@expfun, start_point);
Aest = est(1);
Best = est(2);
Cest = est(3);

Now the optimization arrives at a stable estimate, because the computation is deterministic:

Aest = -0.1266
Best =  1.5106
Cest = 10.2314

The SSE of this estimate is 0.004041.

This is what the data (blue dots) and fitted curves (green: guess, red: optimized) look like:

Answer 2

Here is the whole function in all its glory - special thanks to A. Donda!

function [xval, yval] = curve_fitting_exponential_1_optimized(x,y,xval)

x = x(:);  % bring the data into the standard, more convenient format of column vectors
y = y(:);

Aguess = min(y) - (max(y)-min(y)) / 2;
guess = [ones(size(x)), -x] \ log(y - Aguess);
Bguess = exp(guess(1));
Cguess = guess(2);

start_point = [Aguess, Bguess, Cguess];
est = fminsearch(@expfun, start_point);

    function [sse, FittedCurve] = expfun(params)
        A = params(1);
        B = params(2);
        C = params(3);
        FittedCurve = A + B .* exp(-C * x);
        ErrorVector = FittedCurve - y;
        sse = sum(ErrorVector .^ 2);
    end

yval = est(1)+est(2) * exp(-est(3) * xval);
end