Problem with non-linear least squares fit to a non-linear model function using Gauss-Newtons method

Hi I'm writing a matlab code that will determine least squares fit to a non-linear model function using Gauss-Newtons method. The problem is that what I get in the end is not a good fit and I'm getting a lot of these warnings “Warning: Rank deficient, rank = 3”.

In the figure below the experimental data is plotted as ‘+’ and the red line is the fitted function.

This is how it's supposed to look like.

Below is the function I use to determine partial derivatives of the model function.

function [p1fun,p2fun,p3fun,p4fun] = derivat1(mfun)
syms x a1 a2 a3 a4
p1 = char(diff(mfun(x,a1,a2,a3,a4),a1));
p2 = char(diff(mfun(x,a1,a2,a3,a4),a2));
p3 = char(diff(mfun(x,a1,a2,a3,a4),a3));
p4 = char(diff(mfun(x,a1,a2,a3,a4),a4));
str1 = strcat('@(x,a1,a2,a3,a4)',p1);
p1fun = str2func(str1);
str2 = strcat('@(x,a1,a2,a3,a4)',p2);
p2fun = str2func(str2);
str3 = strcat('@(x,a1,a2,a3,a4)',p3);
p3fun = str2func(str3);
str4 = strcat('@(x,a1,a2,a3,a4)',p4);
p4fun = str2func(str4);
clear x a1 a2 a3 a4
end

And this is the main script for the program.

%----------gnNonlinear-----------
clc,clear
load activity
mfun = @(x,a1,a2,a3,a4) a1*exp(-a2*x)+a3*exp(-a4*x); %function for radioactive activity from a sample containing two radioactive isotopes
[p1fun,p2fun,p3fun,p4fun] = derivat1(mfun);
n = length(t);
%this are the initial guess of the parameters to optimize
a1 = 5000;
a2 = 0.05;
a3 = 2000;
a4 = 0.1;
niter = 1;
tol = 1E-10;
norm = 1;
p1 = zeros(n,1);
p2 = zeros(n,1);
p3 = zeros(n,1);
p4 = zeros(n,1);
while norm>tol && niter<50
    for i = 1:n
        p1(i) = p1fun(t(i),a1,a2,a3,a4);
        p2(i) = p2fun(t(i),a1,a2,a3,a4);
        p3(i) = p3fun(t(i),a1,a2,a3,a4);
        p4(i) = p4fun(t(i),a1,a2,a3,a4);
    end
    A = w.*[p1 p2 p3 p4];
    ytilde = y - mfun(t,a1,a2,a3,a4);
    da = A\(w.*ytilde);
    norm = sqrt(da(1).^2 + da(2).^2 + da(3).^2 + da(4).^2);
    a1 = a1 + da(1);
    a2 = a2 + da(2);
    a3 = a3 + da(3);
    a4 = a4 + da(4);
    niter = niter + 1;
end
disp('a1'),disp(a1)
disp('a2'),disp(a2)
disp('a3'),disp(a3)
disp('a4'),disp(a4)
disp('niter'),disp(niter)
plot(t,y,'+')
hold on
xm = 0:0.1:100;
ym = a1*exp(-a2*xm)+a3*exp(-a4*xm);
plot(xm,ym,'r')
xlabel('x'),ylabel('y')

The file activity.mat has t, y and w vectors. t = time and w = weight.

Why is it that I'm not getting a better fit? Where might I have done wrong?

This is the theory and pseudo code for the method.

Answer 1

[Alan Weiss]

I suggest that you try to follow the example Curve Fitting via Optimization. If you have an Optimization Toolbox license, you would do better to follow the example Nonlinear Data-Fitting.

Alan Weiss

MATLAB mathematical toolbox documentation