I am confused about why my code doesn't lead to a fitting result.

x=[0 1 2 3 4 5 6 7 8 9 10]

x = 1×11

0 1 2 3 4 5 6 7 8 9 10

xs = 0:0.01:120;
y=[0 4 9 16 25 36 49 64 81 100 121]

y = 1×11

0 4 9 16 25 36 49 64 81 100 121

myfunc=@fitting;
rng('shuffle');
best_r = inf;
best_p = [-inf, -inf];
best_c = [-inf];
a=0.00001;
b=3;
d=2;

for iters = 1 : 100
    p0 = a+(b-a).*rand(1,2);
    c0 = a+(d-a).*rand(1,1);
    ytotal0 = fitting(p0,c0,xs);
    try
        [p, c, r] =nlinfit(x,y,myfunc,p0,c0);
        if any(p < 0), or (c < 0)
            fprintf('rejected negative output on iteration %d\n', iters);
            continue;
        end
        
        if any(norm(r)==best_r)
            continue;
        end
        
        nr = norm(r);
        
        
        if nr < best_r
            fprintf('improved on iteration %d\n', iters);
            best_r = nr;
            best_p = p;
            best_c = c;
            
        end
    catch ME
        fprintf('failed iteration %d\n', iters);
    end
end

failed iteration 1
failed iteration 2
failed iteration 3
failed iteration 4
failed iteration 5
failed iteration 6
failed iteration 7
failed iteration 8
failed iteration 9
failed iteration 10
failed iteration 11
failed iteration 12
failed iteration 13
failed iteration 14
failed iteration 15
failed iteration 16
failed iteration 17
failed iteration 18
failed iteration 19
failed iteration 20
failed iteration 21
failed iteration 22
failed iteration 23
failed iteration 24
failed iteration 25
failed iteration 26
failed iteration 27
failed iteration 28
failed iteration 29
failed iteration 30
failed iteration 31
failed iteration 32
failed iteration 33
failed iteration 34
failed iteration 35
failed iteration 36
failed iteration 37
failed iteration 38
failed iteration 39
failed iteration 40
failed iteration 41
failed iteration 42
failed iteration 43
failed iteration 44
failed iteration 45
failed iteration 46
failed iteration 47
failed iteration 48
failed iteration 49
failed iteration 50
failed iteration 51
failed iteration 52
failed iteration 53
failed iteration 54
failed iteration 55
failed iteration 56
failed iteration 57
failed iteration 58
failed iteration 59
failed iteration 60
failed iteration 61
failed iteration 62
failed iteration 63
failed iteration 64
failed iteration 65
failed iteration 66
failed iteration 67
failed iteration 68
failed iteration 69
failed iteration 70
failed iteration 71
failed iteration 72
failed iteration 73
failed iteration 74
failed iteration 75
failed iteration 76
failed iteration 77
failed iteration 78
failed iteration 79
failed iteration 80
failed iteration 81
failed iteration 82
failed iteration 83
failed iteration 84
failed iteration 85
failed iteration 86
failed iteration 87
failed iteration 88
failed iteration 89
failed iteration 90
failed iteration 91
failed iteration 92
failed iteration 93
failed iteration 94
failed iteration 95
failed iteration 96
failed iteration 97
failed iteration 98
failed iteration 99
failed iteration 100

c = best_c

c = -Inf

p = best_p

p = 1×2

-Inf -Inf

%c = lsqcurvefit(myfunc,c0,x,y)
ytotal = fitting(p,c,x);

Warning: Failure at t=0.000000e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (7.905050e-323) at time t.

%for adding points
hold on
plot(x,y,'o', 'displayname', 'original')

hold off

xlim auto
ylim auto
legend show

%for adding points
hold on
xs = 0:0.1:150;
ytotal = fitting(p,c,xs);

Warning: Failure at t=0.000000e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (7.905050e-323) at time t.

plot(xs,ytotal)
hold off

R2=1-(sum((fitting(p,c,x)-y).^2)/sum((y-mean(y)).^2));

Warning: Failure at t=0.000000e+00. Unable to meet integration tolerances without reducing the step size below the smallest value allowed (7.905050e-323) at time t.

function ytotal = fitting(p,c,x)
xspan = x;
y0 = zeros(2,1);
[t,ye] = ode15s(@(x,y)eq2(x,y,p,c), xspan, y0);
%protect against failure of integration
if t(end) ~= xspan(end)
    ytotal = inf(numel(xspan),1);
else
    ytotal = ye(:,1)+ye(:,2)+c(1);
end
end 
function dy=eq2(x,y,p,c)
%y(1)=CE,y(2)=LE
dy=zeros(2,1);
dy(1) = p(1) .* x;
dy(2) = p(2);


end

Answer 1

[Walter Roberson]

I had to fix every different things to get it to work.

The improvement to summarize failed iterations instead of listing each one was cosmetic, not strictly required. It took a while before the iterations started working properly and I was tired of the long stream of fail messages.

x=[0 1 2 3 4 5 6 7 8 9 10]

x = 1×11

0 1 2 3 4 5 6 7 8 9 10

xs = 0:0.01:120;
y=[0 4 9 16 25 36 49 64 81 100 121]

y = 1×11

0 4 9 16 25 36 49 64 81 100 121

rng('shuffle');
best_r = inf;
best_p = [-inf, -inf];
best_c = [-inf];
a=0.00001;
b=3;
d=2;

arefailing = false;
firstfail = -inf;
currentfail = -inf;

for iters = 1 : 100
    p0 = a+(b-a).*rand(1,2);
    c0 = a+(d-a).*rand(1,1);
    ytotal0 = fitting(p0,c0,xs);
    try
        [pc, r] = nlinfit(x,y(:),@(pc,x)fitting(pc(1:2),pc(3),x),[p0,c0]);
        %with the try/catch we never reach this statement if the nlfit
        %succeeded, so the fact we got here means that an iteration worked
        %(but might have returned values we do not want.)
        if arefailing
                fprintf('failed iterations %d to %d\n', firstfail, currentfail);
                arefailing = false;
        end
        p = pc(1:2); c = pc(3);
        if any(p < 0), or (c < 0)
            fprintf('rejected negative output on iteration %d\n', iters);
            continue;
        else
            
        end
        
        if any(norm(r)==best_r)
            continue;
        end
        
        nr = norm(r);
        
        
        if nr < best_r
            fprintf('improved on iteration %d\n', iters);
            best_r = nr;
            best_p = p;
            best_c = c;
            
        end
    catch ME
        if arefailing
            currentfail = iters;
        else
            arefailing = true;
            firstfail = iters;
            currentfail = iters;
            disp(ME)
            if ~isempty(ME.cause); disp(ME.cause{1}); end
        end
    end
end

improved on iteration 1
improved on iteration 4
improved on iteration 19

if arefailing
  fprintf('failed iterations %d to %d\n', firstfail, currentfail);
  if firstfail == 1   %everything failed
      fprintf('failed all iterations!!! Not doing any plotting\n');
  else
      arefailing = false;   %we got at least one valid iteration
  end
end

if ~arefailing
    c = best_c
    p = best_p
    
    
    
    %c = lsqcurvefit(myfunc,c0,x,y)
    ytotal = fitting(p,c,x);
    %for adding points
    hold on
    plot(x,y,'o', 'displayname', 'original')
    
    hold off
    
    xlim auto
    ylim auto
    legend show
    
    %for adding points
    hold on
    xs = 0:0.1:150;
    ytotal = fitting(p,c,xs);
    plot(xs,ytotal)
    hold off
    
    R2=1-(sum((fitting(p,c,x)-y).^2)/sum((y-mean(y)).^2));
end

c = 0.4196

p = 1×2

1.9650 2.2203

function ytotal = fitting(p,c,x)
    xspan = x;
    y0 = zeros(2,1);
    [t,ye] = ode15s(@(x,y)eq2(x,y,p,c), xspan, y0);
    %protect against failure of integration
    if t(end) ~= xspan(end)
        ytotal = inf(numel(xspan),1);
    else
        ytotal = ye(:,1)+ye(:,2)+c;
    end
end 
function dy=eq2(x,y,p,c)
    %y(1)=CE,y(2)=LE
    dy=zeros(2,1);
    dy(1) = p(1) .* x;
    dy(2) = p(2);
end

[Jinglei]

Thank you for your help but the fitting was still a failure with no fitting curve obtained.

[Walter Roberson]

original x are 0 to 10. The plot for them is the blue lower left. The rest of the plot is extrapolation to 120, your xs

[Jinglei]

That's right. And I tried changing the range of xs to (0,10), which, however, led to the same fitting result. In addition, when I tried the same code for more parameters (eight for p and two for c, that is, p=pc(1:8), c=pc(9:10), and of course, for a different model), it didn't come to any result again. In fact, it got stuck in the middle this time, running forever without stopping manually. Do you know what was wrong?Is it caused by to many parameters or any errors in the running process?

[Walter Roberson]

You do not fit against xs, only against x. xs is used only to display the final plot.

[Walter Roberson]

In fact, it got stuck in the middle this time, running forever without stopping manually. Do you know what was wrong?

No. I do not know what is wrong with code that I have not seen.

[Jinglei]

The code can be seen as below.

x = [0
    5
    10
    15
    20
    25
    30
    35
    40
    45
    50
    55
    60
    65
    70
    80
    90
    100
    110
    120
    ];
xs = 0:0.01:120;
y = [0
    1.415699757
    1.044160175
    1.342848859
    1.357419038
    1.138866343
    1.386559398
    1.452125206
    1.430269937
    1.153436523
    1.26271287
    1.408414667
    1.357419038
    1.430269937
    1.095155804
    1.488550656
    1.503120835
    1.386559398
    1.401129577
    1.532261195
    ];
rng('shuffle');
best_r = inf;
best_p = [-inf, -inf,-inf,-inf,-inf,-inf,-inf,-inf];
best_c = [-inf,-inf];
a=0.00001;
b=1;
d=2;

arefailing = false;
firstfail = -inf;
currentfail = -inf;
for iters = 1 : 100
    p0 = a+(b-a).*rand(1,8);
    c0 = a+(d-a).*rand(1,2);
    ytotal0 = fitting(p0,c0,xs);
    try
        [pc, r] = nlinfit(x,y(:),@(pc,x)fitting(pc(1:8),pc(9:10),x),[p0,c0]);
        %with the try/catch we never reach this statement if the nlfit
        %succeeded, so the fact we got here means that an iteration worked
        %(but might have returned values we do not want.)
        if arefailing
            fprintf('failed iterations %d to %d\n', firstfail, currentfail);
            arefailing = false;
        end
        p = pc(1:8); c = pc(9:10);
        if any(p < 0), or (c < 0)
            fprintf('rejected negative output on iteration %d\n', iters);
            continue;
        else
            
        end
        
        if any(norm(r)==best_r)
            continue;
        end
        
        nr = norm(r);
        
        
        if nr < best_r
            fprintf('improved on iteration %d\n', iters);
            best_r = nr;
            best_p = p;
            best_c = c;
            
        end
    catch ME
        if arefailing
            currentfail = iters;
        else
            arefailing = true;
            firstfail = iters;
            currentfail = iters;
            disp(ME)
            if ~isempty(ME.cause); disp(ME.cause{1}); end
        end
    end
end

if arefailing
    fprintf('failed iterations %d to %d\n', firstfail, currentfail);
    if firstfail == 1   %everything failed
        fprintf('failed all iterations!!! Not doing any plotting\n');
    else
        arefailing = false;   %we got at least one valid iteration
    end
end
if ~arefailing
    c = best_c
    p = best_p
    
    
    
    %c = lsqcurvefit(myfunc,c0,x,y)
    ytotal = fitting(p,c,x);
    %for adding points
    hold on
    plot(x,y,'o', 'displayname', 'original')
    
    hold off
    
    xlim auto
    ylim auto
    legend show
    
    %for adding points
    hold on
    xs = 0:0.1:150;
    ytotal = fitting(p,c,xs);
    plot(xs,ytotal)
    hold off
    
    R2=1-(sum((fitting(p,c,x)-y).^2)/sum((y-mean(y)).^2));
end

function ytotal = fitting(p,c,x)
xspan = x;
y0 = zeros(4,1);
[t,ye] = ode15s(@(x,y)eq2(x,y,p,c), xspan, y0);
%protect against failure of integration
if t(end) ~= xspan(end)
    ytotal = inf(numel(xspan),1);
else
    ytotal = ye(:,1)+ye(:,2)+ye(:,3)+ye(:,4);
end
end 
function dy=eq2(x,y,p,c)
%y(1)=CE,y(2)=LE
dy=zeros(4,1);
dy(1) = p(1) .* (2.5 - y(1)-y(2)-y(3)-y(4)) .* (c(1)-y(1)-y(2)-y(4)) - p(2) .* y(1);
dy(2) = p(3) .* y(1)-p(4).*y(2);
dy(3) = p(5).*(2.5-y(1)-y(2)-y(3)-y(4)).*(c(2)-y(3)-y(4))-p(6).*y(3);
dy(4) = p(7).*y(3).*(c(2)-y(1)-y(2)-y(4))-p(8).*y(4);

end

Answer 2

[Bora Eryilmaz]

If you print out the exception inside the catch block, it will tell you what the error is.

 
catch ME
     ME
     fprintf('failed iteration %d\n', iters);
     

To start with, you are not passing the right number of arguments to the fitting() function: it would be passed two arguments from nlinfit, but your code expects three input arguments.

[Jinglei]

Thank you for your reply.

Firstly, I think those warnings were caused by the range of parameters I used, in which I gave the initial values as inf. I tried exact values instead and the warnings didn't appear. However, this doesn't help for the fitting. At the same time, I have found that the final values for all parameters are the same as the initially given ones (c_best, p_best, and r_best). In other words, those 'if' codes are not runned, which is the real problem.

In addition, as for the number of arguments you mentioned, I am not sure which three arguments you were talking about. If the third one is r, then I don't think it matters since it is not related to the fitting process but just the deviation to get the optimal simulation.

[Bora Eryilmaz]

You function fitting() takes three input arguments, p, c, and x. However, I think nlinfit will call that function with only two arguments.

Also, please read the documentation for nlinfit: https://www.mathworks.com/help/stats/nlinfit.html. Among other things, nlinfit does not return the estimated coefficients individually; it returns them together as a vector.

Finally, using Inf as initial guesses of parameters does not sound right; initial parameters are supposed to be at reasonable, ballpark values.

[Jinglei]

Thank you for your reply.

Firstly, I think what you meant is that the nlinfit can only simulate two arguments including x, that is, only one parameter could be included, c or p. However, this doesn't work for my situation, even if I combined c and p together to be one, still, no fitting result could be achieved.

In addition, with the 'if' codes not running, c and p will be the initial values given instead of that obtained after the for codes. This means running nlinfit correctly or not would not matter since p and c will be assigned as p_best and c_best as intially given. Therefore, I would like to know what makes 'if' codes not working.

Further, it is strange that each run of the "for" code gave falure (100 falures as can be seen) in the runing process and I would like to know how to fix that.

[Walter Roberson]

For the case where you are passing in fixed values (for example a table of oxygen absorption) then see http://www.mathworks.com/help/matlab/math/parameterizing-functions.html for how to construct a function handle that the function would be happy with.

The objective function you pass to nlinfit must expect exactly two parameters: the vector of proposed coefficients, and the set of X values that the proposed parameters are to be applied to. If you have a model in which you are wanting more than one model variable to be altered, then put all of the model variables into a vector. For example if you have a model in parameters c and p and you want nlinfit to alter both of those, then create

function y = fitting(cp, x)
    c = cp(1);
    p = cp(2);
    y = appropriate expression
end

[Jinglei]

Thank you for your suggestion. However, in this case, c and p are taken as un-defined parameters by Matlab. And I am confused how to fix that. Could you show me that?

[Jinglei]

Thank you for your help and it does work this time. The only problem is the obtained parameters to some extent away from the optimal ones, which may not be a big problem for my data. Thank you for your patience.