lsqcurvefit residuals don't agree with the manual calculation

Hello,

I'm using lsqcurvefit to fit experimental modeling data to the analytic solution of a solute transport equation. Unfortunately, the calculation of the residuals via lsqcurvefit, do not match my manual calculation of the residuals between the best-fit curve and experimental modeling data.

I call lsqcurvefit like this:

dt = table2array(outletConcentrations(:,1)); % first column of this table is a time vector
allconcs = [concA, concB, concC]; concAvg = mean(allconcs,2); % these columns are concentration vectors
xmids = [0.01999, 0.01999, 0.01999]; % [m] mid point fluid sampling format [centerA, CenterB, CenterC]
hardparams = [xmids(1), V_out2p2p, 480]; % these are intentionally hardwired parameters 
initialparams = [2e-1, 0.1/0.5]; % alpha, beta -- these are the parameters I'd like to optimize using lsqcurvefit


[params_est, resnorm,residual,exitflag,output] =...
    lsqcurvefit(@(initialparams,dt) soln1dNoDisp(initialparams,hardparams,dt),initialparams,dt,concAvg);
% calling the function sol1dnodisp as such, produces the output from sol1dnodisp which has the following syntax
% function [c,cim,alpha,beta]=soln1dNoDisp(initialparams,hardparams,dt)
 
% Also, note that c is the array being fit to concAvg by adjusting alpha
% and beta. 

I have lsqcurvefit outputting the residuals for error checking, however when integrating the vector of residuals I get pretty large numbers (-135 mol/m^3). So I decided to check the work with a manual calculation of the residuals and their displacement as follows:

LDresiduals = c - concAvg; % The help menu of lsqcurvefit indicates the residuals are calculated as: 
% FUN(X,XDATA)-YDATA. I believe this is exactly what I have done manually here.
displaced = LDresiduals - mean(LDresiduals-residual);

You can see in the attached image that the residuals seem to be shifted in magnitude. The 'residuals, displaced' data is my checking to see if they are displaced, which seems likely true within the error produced by machine precision (mean(displaced-residual) = 1.99e-17). Is anyone able to explain why I'm seeing this?

I've attached the experimental model data and best fit curve data for reference; 'c mobile' is being fit to 'average outlet'

[Torsten]

We don't see how you compute "c" from "params_est".

Further make sure that you return "c" in "soln1dNoDisp", not "c - concAvg".

Your calculated "LDresiduals" should equal "residual" returned by lsqcurvefit.

[Charles Dorchester]

c is calculated in the function sol1dnodisp, the definition of which was given in a comment. Here it is again:

function [c,cim,alpha,beta]=soln1dNoDisp(initialparams,hardparams,dt).

Unfortunately I'm not able to share the sol1dNoDisp code. However I have provided figures of the output. See the bottom photo which shows 'c', 'cim', and 'concAvg'.

Yes, I see that 'LDresiduals' and 'residual' are not equal, but I'm aware that they should be equal.

I think I'm probably just confused on how lsqcurvefit should be called.

if I'm working with concentration breakthrough curves with respect to time, my output should be f(c,t). Therefore, if I'm following the example given by matlab:

% Matlab's example:
xdata = [3; 1; 4];           % example xdata
ydata = 6*exp(-1.5*xdata)+3; % example ydata     
c = -1.5;                    % define parameter 
x = lsqcurvefit(@(x,xdata) myfun(x,xdata,c),[5;1],xdata,ydata)

Does that mean the following variables should translate from the example to my problem? :

xdata -> time

x -> concentration

ydata -> concentration from the best fit curve

[Torsten]

Does that mean the following variables should translate from the example to my problem? :

xdata -> time

x -> concentration

ydata -> concentration from the best fit curve

No.

xdata -> time

x -> parameters supplied from lsqcurvefit to myfun which are to be used to calculate the simulated concentrations by using ode15s

ydata -> concentrations from measurements which are to be reproduced by your model

Answer 1

[Charles Dorchester]

Thanks for the help.

I was able to sort it out when I realized that the output of my function was producing results between 0 and 1. I had then scaled these results to fit the experimental modeling data within the script, but outside of the function.

Scaling the results within the function produces the correct results when it is called by lsqcurvefit.