[Help-gsl] strange NLS behaviour

[viktor drobot]

Hi, community!

I model kinetics of enzymatic reactions. There are problems for which it's
necessary to find suitable values of kinetic parameters so the real kinetic
data could be described by our model.

To facilitate my test runs I generate some "real kinetic data" by solving
the system of ODEs with known kinetic parameters and then add some Gaussian
noise. This is the reference curves. After that I try to solve the inverse
problem - find out the most appropriate values of these parameters with
initial approximation that is close enough to the true values. I use GLS
NLS routines (v 2.2.1) for that.

First of all I solve ODEs with another approximation and then find the
differences (Y[t] - y[t]) between reference points Y[t] (generated as
stated above) and the new calculated y[t], where 't' is the time point.
Usually I process several kinetic curves at once. I store these differences
in 'f' vector of function to be minimized. Because of the nature of the
problem there is no analytical Jacobian possible so I set corresponding
field of gsl_multifit_nlinear_fdf structure to NULL. In case of unweighted
NLS I multiply the diagonal elements of the covariance matrix by the
variance of the residuals. Then the NLS solver does all work for me and I
get some strange results.

The output says that there was no Jacobian evaluations. Kinetic parameters
obtained seems to be quiet right but other statistics weird enough. Look at
the example below.

InverseSolver v0.1.0
********************
Run ended at Thu Sep  8 17:13:27 2016


====================
Optimized parameters
--------------------
Kcat = 2.4849549376350680e+01 +/- 2.3301006770134422e-01 (true value is
2.5e+01)
Km   = 2.5962066490476204e-05 +/- 1.2805709675629306e-06 (true value is
2.5e-05)
Kp   = 5.6728528418318680e-05 +/- 3.3698382714302584e-06 (true value is
5.0e-05)


====================
Minimization details
--------------------
The fit converged after 13 iteration(s)
Reason for stopping: small step size
Objective function evaluations = 69
Number of data points          = 380
Model parameters               = 3
dof                            = 377
chi                            = 2.7339109755226917e-05
chisq                          = 7.4742692220834363e-10
chisq/dof                      = 1.9825647803934845e-12
Time elapsed: 0.05957 secs


====================================
Correlation matrix of fit parameters
------------------------------------
     Kcat        Km          Kp
Kcat  1.0000e+00  8.5928e-01  6.5048e-01
Km    8.5928e-01  1.0000e+00  9.4124e-01
Kp    6.5048e-01  9.4124e-01  1.0000e+00


==================
Iterations history
------------------
Iter # Kcat                    Km
Kp                      cond(J)                 chi
0       2.0000000000000000e+01  2.0000000000000002e-05
2.0000000000000002e-05                     inf  3.1567955122657597e-04
1       2.7173704932869438e+01  2.6087621471921223e-05
3.0462571586587205e-05  2.9025496731552132e+06  8.3486976629633576e-05
2       2.4543231868632770e+01  2.5481466984423840e-05
4.6317587677535904e-05  3.0659026221618010e+06  4.9189276130628142e-05
3       2.4799130142472123e+01  2.5333806359622357e-05
5.3688047773966246e-05  2.1896626012014118e+06  2.7592054291596016e-05
4       2.4830032619351048e+01  2.5841007936620306e-05
5.6370153753761361e-05  2.0525957706617371e+06  2.7340238653871502e-05
5       2.4848466472445963e+01  2.5956371903529538e-05
5.6716387131234307e-05  2.0040307592909317e+06  2.7339110581465998e-05
6       2.4849535469371190e+01  2.5961998929997054e-05
5.6728382061723916e-05  1.9990772957287871e+06  2.7339109755370303e-05
7       2.4849547296667978e+01  2.5962054330720485e-05
5.6728499949112989e-05  1.9989372038075817e+06  2.7339109755230966e-05
8       2.4849547539278458e+01  2.5962056427776209e-05
5.6728505891823464e-05  1.9989391076721402e+06  2.7339109755230180e-05
9       2.4849549814327389e+01  2.5962073280806688e-05
5.6728551401861829e-05  1.9989361463323494e+06  2.7339109755228415e-05
10      2.4849550175380344e+01  2.5962072172440512e-05
5.6728543239809625e-05  1.9989390874822082e+06  2.7339109755227405e-05
11      2.4849549199258245e+01  2.5962065484919628e-05
5.6728525792610483e-05  1.9989382142781706e+06  2.7339109755227290e-05
12      2.4849549376350691e+01  2.5962066490476238e-05
5.6728528418318701e-05  1.9989384807737938e+06  2.7339109755226951e-05
13      2.4849549376350680e+01  2.5962066490476204e-05
5.6728528418318680e-05  1.9989391375842493e+06  2.7339109755226917e-05


=============================
Per-file Fischer's statistics
-----------------------------
[/home/viktor/science/enzyme_kinetics/data/test/172-186/0172_1s-0075.dat]
                  dof = 51
Number of data points = 54
         F-statistics = 3.2567463900742809e+04

[/home/viktor/science/enzyme_kinetics/data/test/172-186/0173_1s-0100.dat]
                  dof = 62
Number of data points = 65
         F-statistics = 1.0932769130444922e+04

[/home/viktor/science/enzyme_kinetics/data/test/172-186/0174_1s-0125.dat]
                  dof = 72
Number of data points = 75
         F-statistics = 2.2345670716799857e+04

[/home/viktor/science/enzyme_kinetics/data/test/172-186/0175_1s-0150.dat]
                  dof = 84
Number of data points = 87
         F-statistics = 4.4888120828034567e+04

[/home/viktor/science/enzyme_kinetics/data/test/172-186/0176_1s-0175.dat]
                  dof = 96
Number of data points = 99
         F-statistics = 3.9235545459615321e+04

chisq/dof value is too small. It's much smaller than 1. Condition number of
Jacobian is too big, much bigger than 1. Error estimates are too small
despite of some Gaussian noise in my generated data. Seems like the
calculation is unstable.

I completely stuck with these strange results. What I'm doing wrong? Any
help is greatly appreciated!

-- 
С уважением,
Дробот Виктор

plot_c.png
Description: PNG image