how to modify weights of each square error in the SSE objective function of lsqcurvefit ?

The lsqcurvefit uses the sum of square of errors (SSE) minimization, right? I have tried using it using a set of chemical reactions where some species concentrations are in the 1E+01 magnitude and some are in 1E-02. What I have found is that the local minimum found by the solver does not give a good fit for the concentrations of species that is in the 1E-2 magnitude. Is it because its square of error does not contribute much to the SSE minimization? If so, how do I modify the objective function to get a local minimum that also fits well with the concentrations that are of the 1E-02 magnitude?

I have followed the https://www.mathworks.com/matlabcentral/answers/43439-monod-kinetics-and-curve-fitting#comment_89455 solution.

Answer 1

[Apeksha Bagrecha]

Here is an example to depict how to use weights, explaining how using the "lsqnonlin" function instead of the "lsqcurvefit" will help you achieve this.

I am attaching two files with this email: "fun.m" and "sample.m" (which is the main function). Please refer to this code and read the explanation below.

"lsqcurvefit" itself applies to a special case of "lsqnonlin" when the residual vector has the form

 r(x)= F(x,xdata)-ydata

Conversely, lsqnonlin can be applied to any (differentiable) residual function r(x).

In our case, since we want to apply a "weights" vector, we'll define weighted versions of the data and the model function, then use nonlinear least squares to make the fit:

function F = fun(xdata, ydata, w, ycurve, a)

ycurvedata = ycurve(a);

F = sqrt(w.*(ycurvedata - ydata).^2);

end

Here, "ycurve" refers to the function handle for the function that you are trying to fit, while "w" is the weights vector. The above code snippet is contained in the function file "fun.m".

We now solve this nonlinear least-squares (nonlinear data-fitting) problem using lsqnonlin:

xdata = [1 2 3 5 7 10]';

ydata = [109 149 149 191 213 224]';

w = [1, 1, 5, 5, 5, 5]';

ycurve = @(a) a(1).*(1-exp(-a(2).*xdata));

x0 = [240; 0.5];

F = @(a) fun(xdata, ydata, w, ycurve, a);

[x,resnorm,residual,exitflag,output,lambda,jacobian] = lsqnonlin(F, x0);

The above code snippet is present in "sample.m". Subsequently, we compute the fitted response values confidence intervals, and plot it.

When you run "sample.m", you will notice that the two down-weighted points are not fit as well by the curve as the remaining points. That is as you would expect for a weighted fit.

You can adapt this example for your own workflow.

[Arniel D]

Wow! This is a good way to do it. Thank you so much! Btw, where can I access the attached file you mentioned? I don't see it.

[Apeksha Bagrecha]

I have added the files

[Arniel D]

Hi, I'm having difficulty with implimenting this to my system. The weights applied here are for each xdata point and there's only one function, y=f(x,a). I am dealing with many species with their own nonlinear concentration profiles (C vs t).

My problem is similar to the Monod kinetics and curve fitting problem where I have concentration profiles (concentration vs time) for several species and I wanted to put weights on each kinetic rate of species (defined by ODEs) since I found that when I apply the lsqcurvefit, the concentration profile of the intermediate species are not fitted well because their concentrations (y) are much smaller in magnitude with the rest (like reactants and products).

So I have a fixed time, say t = [0, 1, 2, 3, 4, 5, 6, 7, 8]; and I have several species like A, B, C, D, E, F, where each species have their own rate equation (dCi/dt) that includes the parameters I wanted to determine. Each species have experimental data at each time. For example, A = [10, 8, 4, 3, 2, 1, 0.5, 0.4, 0.4]; and so on.

How can I put weights for each species?