lmcurvefit

curve fitting using Levenberg Marquardt algorithm

% the following examples are available here
% https://www.mathworks.com/help/optim/ug/lsqcurvefit.html
% you can compare to the result from lmcurvefit to that of inbuilt matlab
% function lsqcurvefit

%% Example 1 (Unconstrained Curve Fitting)
close all
xdata = [0.9 1.5 13.8 19.8 24.1 28.2 35.2 60.3 74.6 81.3]';
ydata = [455.2 428.6 124.1 67.3 43.2 28.1 13.1 -0.4 -1.3 -1.5]';
x0 = [100;-1];
times = linspace(xdata(1),xdata(end))';
scatter(xdata, ydata, 'o'); hold on;
plt = plot(times, myfun1(x0, times), 'r');
model = @(x,xdata) myfun1(x, xdata, times, plt);
[x_lm, ~, ~, ~, output] = lmcurvefit(model, x0, xdata, ydata,[],[],[],[])

%% Example 2 (Box Constrained Curve Fitting)
close all
xdata = linspace(0, 3)';
ydata = exp(-1.3*xdata)+0.05*rand(size(xdata));
lb = [0;-2];
ub = [3/4; -1];
x0 = [1/2;-2];
scatter(xdata, ydata, 'o'); hold on;
plt = plot(xdata, myfun2(x0, xdata), 'r');
model = @(x,xdata) myfun2(x, xdata, xdata, plt);
[x_lm, ~, resnorm_lm, residual_lm, output] = ...
            lmcurvefit(model, x0, xdata, ydata, [], [], lb, ub);

%% Example 3 (Linear InEquality Constraint)
close all; clc
rng default
xdata = linspace(2,7)';
ydata = myfun3([2,4,5,0.5]',xdata) + 0.1*randn(size(xdata));
lb = zeros(4,1);
ub = 7*ones(4,1);
A = [-1 -1 1 1];
b = 0;
startpt = [1 2 3 1]';
options = optimoptions(@lsqcurvefit, Display='iter');
scatter(xdata, ydata, 'o'); hold on;
plt = plot(xdata, myfun3(startpt,xdata), 'r');
fineq = @(x)A*x - b; 
fun = @(x, xdat) myfun3(x,xdat, plt);
[x_lm, ~, resnorm_lm, residual_lm, output_lm] = ...
    lmcurvefit(fun, startpt, xdata, ydata, fineq, [], lb, ub);

%% Example 4 (Nonlinear InEquality Constraint)
close all; clc
rng default
xdata = linspace(2,7)';
ydata = myfun3([2,4,5,0.5]',xdata) + 0.1*randn(size(xdata));
lb = zeros(4,1);
ub = 7*ones(4,1);
startpt = [1 2 3 1]';
options = optimoptions(@lsqcurvefit, Display='iter');
scatter(xdata, ydata, 'o'); hold on;
plt = plot(xdata, myfun3(startpt,xdata), 'r');
fineq = @(x)x(1)^2 + x(2)^2 - 4^2; 
fun = @(x, xdat) myfun3(x,xdat, plt);
[x_lm, ~, resnorm_lm, residual_lm, output_lm] = ...
    lmcurvefit(fun, startpt, xdata, ydata, fineq, [], lb, ub);

%% model functions
function F = myfun1(x,xdata, times, plt)
    F = x(1)*exp(x(2)*xdata);
    if(nargin > 2)
        plt.YData = x(1)*exp(x(2)*times); 
        drawnow; pause(0.01);
    end
end

function F = myfun2(x,xdata, times, plt)
    F = x(1)*exp(x(2)*xdata);
    if(nargin >2)
        plt.YData = x(1)*exp(x(2)*times); 
        drawnow; pause(0.01);
    end
end

function F = myfun3(x,xdata, plt)
    a = x(1); b = x(2); t0 = x(3); c = x(4);
    F = a + b*atan(xdata - t0) + c*xdata;
    if(nargin>2)
        plt.YData = F; 
        drawnow; pause(0.01);
    end
end