Nonlinear curve fitting issue.

Hello guys,

I am trying to implement a routine for fitting electrophoretic data from my experiments.

The aim is to derive kinetic parameters for the interaction of biomoecules from the relative areas of peaks in the electropherogram, based on the areas of the peaks in the dataset.

Since all relevant differential equations are known and since the set of equations has an analytical solution, as described here:

http://www.yorku.ca/skrylov/Publications/Anal_Chem_2004_v76_pp1507-1512.pdf

I set about entering the relevant equations (6, 8, 13 from the referenced manuscript) in matlab.

The thus created function works and I can use it to simulate electropherograms of interacting species.

Obviuously, I now would like to use the function to fit experimental data and retrieve the parameters (8 in total, Va, Vc, MUa, MUc, k, A0, C0, baseline noise).

Some of these will obviously be correlated.

My problem is, if I supply experimental data and try something like lsqcurvefit:

[x,resnorm,residual] = lsqcurvefit(@(param,xdata) Electropherogram2(param,xdata,column), params0, time, M(1:2400,2),lb, ub, options);

I often get very poor results because I either run out of iterations, I hit some (obviously poorly fitting) local minimum.

Only if I tinker a lot with the starting values and the allowed intervals (i.e. because I know likely values through other experiments) do I end up with more or less decent fits, but teven then, fits are not as good as reported in the original manuscript (fig. 3).

The authors of that manuscript used Excel solver and were kind enough to provide the data used in Fig. 3 but still I cannot seem to end up with fits as good as theirs without nearly literally supplying the nearly correct starting values.

I am not experienced enough to know what I could tweak to make this process less trial-and-error.

Would something like the global optimization toolbox help me?

Any tips are welcome...

Answer 1

[Star Strider]

First, you can deal with the problem of ‘running out of iterations’ by increasing 'MaxIter' and 'MaxFunEvals' in optimset, increasing perhaps both to 1E+5. You might also experiment with decreasing 'TolFun' and 'TolX' to 1E-10 or so. It takes longer, but increases your probability of getting a good fit. You can speed that process by providing an analytic Jacobian, also easy if you have the Symbolic Math Toolbox.

Second, I've found that the trust-region-reflective algorithm lsqcurvefit uses can occasionally find several local minima before it finds the global minimum. I suggest you experiment with 'Algorithm','levenberg-marquardt' in optimset to see if that avoids local minima. The only drawback is that using the Levenberg-Marquardt algorithm doesn't allow you to constrain the parameters. Also, if you have the Statistics Toolbox, it's always a good idea to check your parameters for statistical significance with nlparci. (There is a workaround if you don't have it.)

Third, patternsearch in the Global Optimization Toolbox can be very helpful in situations where an objective function has many local minima. I definitely suggest it, especially if the Levenberg-Marquardt option for lsqcurvefit still gives you problems.

Fourth, objective functions can be difficult to fit if the data vary significantly in magnitude (on the order of 1E±6 or so). The Excel Solver scales automatically (not something I would recommend), but scaling can be helpful in some situations if the parameters affected by scaling are linearly related to the scaled data.

Fifth, nonlinear parameter estimation with complicated problems is invariably trial-and-error. The patternsearch results can depict this vividly.

Sixth, the paper was interesting. I learned something, even though I didn't take the time to read it exhaustively. (I saved it to read later.)

Last, some don't have a high opinion of Excel Solver compared to MATLAB.