Message Boards Message Boards

Global data fitting for a chemical kinetics problem?

Dear all The problem of global data analysis is a common problem among experimenters, especially in chemical kinetics. I have different kinetics (variation of absorbance as a function of time) at different wavelengths. I would like to adjust a set of 10 kinetics. These kinetics are governed by the same kinetic law, the same rate constants (here k3,k4,k5 ) but the absorption coefficients (ea, eh and ef) depend on the wavelength. I usually do these fits with Igor (Wavemetrics) but without any analytical solution I turned to Mathematica. I tried the MultiNonlinearModelFit function 

Clear[k3, k4, k5, ea1, eh1, ef1, ea2, eh2, ef2]
sol = ParametricNDSolve[{Derivative[1][ag][t] == -k3*ag[t], 
    Derivative[1][gox][t] == k3*ag[t], 
    Derivative[1][h][t] == 
     k3*ag[t] - 2*k4*h[t]*h[t] + 2*k5*ff[t]*ff[t], 
    Derivative[1][ff][t] == 2*k4*h[t]*h[t] - 2*k5*ff[t]*ff[t], 
    ag[0] == 0.000006, h[0] == 0, gox[0] == 0, ff[0] == 0}, {ag, gox, 
    h, ff}, {t, 0, 2000}, {k3, k4, k5}];
{agFunc, goxFunc, hFunc, fFunc} = {ag, gox, h, f} /. sol;
model1[k3_, k4_, k5_, ea1_, eh1_, ef1_][
   t_] := (((agFunc[k3, k4, k5][t])*ea1) + (eh1*
      hFunc[k3, k4, k5][t]) + (ef1*fFunc[k3, k4, k5][t]));
model2[k3_, k4_, k5_, ea2_, eh2_, ef2_][
   t_] := (((agFunc[k3, k4, k5][t])*ea2) + (eh2*
      hFunc[k3, k4, k5][t]) + (ef2*fFunc[k3, k4, k5][t]));
fit = MultiNonlinearModelFit[{dat1, 
    dat2}, {model1[k3, k4, k5, ea1, eh1, ef1][t], 
    model2[k3, k4, k5, ea2, eh2, ef2][t]}, {k3, k4, k5, ea1, eh1, ef1, 
    ea2, eh2, ef2}, {t}];
fit["BestFitParameters"]

I have a concern with the output I don't understand. 

ParametricNDSolve::ndsz: At t == 3.127085338387159`, step size is effectively zero; singularity or stiff system suspected.

InterpolatingFunction::dmval: Input value {3.2} lies outside the range of data in the interpolating function. Extrapolation will be used.

...

NonlinearModelFit::nrlnum: The function value {-0.039527+1. f[1.,1.,1.][0.2],-0.0396433+1. f[1.,1.,1.][0.4],-0.0396255+1. f[1.,1.,1.][0.6],-0.0394311+1. f[1.,1.,1.][0.8],-0.0393989+1. f[1.,1.,1.][1.],-0.0394361+1. f[1.,1.,1.][1.2],-0.0394022+1. f[1.,1.,1.][1.4],-0.039186+1. f[1.,1.,1.][1.6],-0.0392537+1. f[1.,1.,1.][1.8],-0.0392642+1. f[1.,1.,1.][2.],<<32>>,-0.0384155+1. f[1.,1.,1.][8.6],-0.0383621+1. f[1.,1.,1.][8.8],-0.0385297+1. f[1.,1.,1.][9.],-0.0384413+1. f[1.,1.,1.][9.2],-0.0384188+1. f[1.,1.,1.][9.4],-0.0382994+1. f[1.,1.,1.][9.6],-0.0383953+1. f[1.,1.,1.][9.8],-0.0381658+1. f[1.,1.,1.][10.],<<850>>} is not a list of real numbers with dimensions {900} at {k3,k4,k5,ea1,eh1,ef1,ea2,eh2,ef2} = {1.,1.,1.,1.,1.,1.,1.,1.,1.}.

Yet it provides me with parameters which for some are far from the expected values: {k3 -> 76620.6, k4 -> -23293.5, k5 -> 10237.6, ea1 -> 169.232, eh1 -> 81.862, ef1 -> 14461.1, ea2 -> 157.347, eh2 -> 80.0374, ef2 -> 14761.5}

Parts of the data are 

dat2 = {{0.2, 0.0366929}, {0.4, 0.0368321}, {0.6, 0.0368786}, {0.8, 
   0.0367551}, {1, 0.036677}, {1.2, 0.036687}, {1.4, 0.0365921}, {1.6,
    0.0365438}, {1.8, 0.0364348}, {2, 0.0364835}, {2.2, 
   0.036489}, {2.4, 0.0363618}, {2.6, 0.0361545}, {2.8, 
   0.0362492}, {3, 0.0362564}, {3.2, 0.0362539}, {3.4, 
   0.0361269}, {3.6, 0.0362652}, {3.8, 0.0361286}, {4, 
   0.0361022}, {4.2, 0.0360444}, {4.4, 0.0360916}, {4.6, 
   0.0360261}, {4.8, 0.0360464}, {5, 0.0359304}, {5.2, 
   0.0358969}, {5.4, 0.0359568}, {5.6, 0.036002}, {5.8, 
   0.0358563}, {6, 0.0360431}, {6.2, 0.0358403}, {6.4, 
   0.0358596}, {6.6, 0.0358963}, {6.8, 0.0359946}, {7, 
   0.0358273}, {7.2, 0.0359167}, {7.4, 0.0357477}, {7.6, 
   0.0357798}, {7.8, 0.0358253}, {8, 0.0358452}, {8.2, 
   0.0358645}, {8.4, 0.035836}, {8.6, 0.0358447}, {8.8, 
   0.0357238}, {9, 0.0358621}, {9.2, 0.0358738}, {9.4, 
   0.0357445}, {9.6, 0.0356904}, {9.8, 0.0357929}, {10, 0.0355983}}

and 

dat1 = {{0.2, 0.039533}, {0.4, 0.0396493}, {0.6, 0.0396315}, {0.8, 
   0.0394371}, {1, 0.0394049}, {1.2, 0.0394421}, {1.4, 
   0.0394082}, {1.6, 0.039192}, {1.8, 0.0392597}, {2, 
   0.0392702}, {2.2, 0.039207}, {2.4, 0.0391561}, {2.6, 
   0.0390186}, {2.8, 0.038951}, {3, 0.0390279}, {3.2, 
   0.0389623}, {3.4, 0.0389681}, {3.6, 0.0389333}, {3.8, 
   0.0387812}, {4, 0.0388705}, {4.2, 0.0388193}, {4.4, 
   0.0387575}, {4.6, 0.0386211}, {4.8, 0.0388169}, {5, 
   0.0387259}, {5.2, 0.038528}, {5.4, 0.0387643}, {5.6, 
   0.03884}, {5.8, 0.0385277}, {6, 0.0385109}, {6.2, 0.0386881}, {6.4,
    0.0385117}, {6.6, 0.038618}, {6.8, 0.0386693}, {7, 
   0.0385695}, {7.2, 0.0385362}, {7.4, 0.0384269}, {7.6, 
   0.038459}, {7.8, 0.038489}, {8, 0.0384817}, {8.2, 0.0384481}, {8.4,
    0.0384566}, {8.6, 0.0384215}, {8.8, 0.0383681}, {9, 
   0.0385357}, {9.2, 0.0384473}, {9.4, 0.0384248}, {9.6, 
   0.0383054}, {9.8, 0.0384013}, {10, 0.0381718}}

other question, but perhaps not far from the problem: is it possible to reduce the search area for the right parameters by specifying an existence interval for each parameter. E.g. 0.001<k3<0.01, 1<k4<400, 500<ea, etc?

Thank you in advance for your advice and help.

Hélène
POSTED BY: Hélène PASQUIER
5 Replies
Sort By:
Jim Baldwin
Jim Baldwin, Retired
Posted 4 years ago

I think the main issue is that both models are overparameterized for the available data (and secondarily better starting values are needed). A plot of the data shows just slight curvature over values of $t$ for both datasets suggesting that maybe just 3 parameters for each curve might be warranted. Or at least that having 6 parameters per curve is too much for the available data.

I want to emphasize that I'm not saying the models are wrong. It's just that the available data (meaning values of $t$ from 0 to 10) doesn't allow for the estimation of all of the parameters. Maybe if there are more than two datasets, then more parameters might be able to be estimated.

Here is what works:

fit = MultiNonlinearModelFit[{dat1, dat2}, {model1[k3, k4, k5, ea1, eh1, ef1][t], 
    model2[k3, k4, k5, ea2, eh2, ef2][t]},
   {{k3, 0.22}, {k4, 0}, {k5, 0}, {ea1, 6600}, {eh1, 6400}, {ef1, 0}, {ea2, 6200}, {eh2, 6000}, {ef2, 0}}, {t}];
mle = fit["BestFitParameters"]
(* {k3 -> 0.2188, k4 -> -5.24086, k5 -> 0., ea1 -> 6623.65, 
 eh1 -> 6364.94, ef1 -> 0.000336612, ea2 -> 6152.56, eh2 -> 5927.64, ef2 -> 0.000419775} *)
fit["CorrelationMatrix"] // MatrixForm

Correlation matrix

fit["ParameterTable"]

Parameter table

Note that the P-values suggest that k4, k5, ef1, and ef2 are not statistically significantly different from zero. However, the huge standard errors suggest that those parameters are not estimated very well. In short, there's not much one can say about those 4 parameters.

Show[ListPlot[{dat1, dat2}],
 Plot[{model1[k3, 0, 0, ea1, eh1, 0][t] /. mle,
   model2[k3, 0, 0, ea2, eh2, 0][t] /. mle}, {t, 0, 10}]]

Plot of data and fit

While MultiNonlinearModelFit is very handy for multiple datasets with common parameters, it does assume that the variability about the individual curves is the same (although with the Weight option one could allow for knowing that the variability is, say, twice as much in one dataset as another but I suspect that is a rare occurrence.) I think an explicit statement in the documentation about the assumption of the equality of the variances is warranted.
POSTED BY: Jim Baldwin

Indeed, you're right, I only sent you the beginning of the data : the data goes up to 850 and not 10.

I managed to do again what you did. Thank you so much for that breakthrough.I guess {k3, 0.22} means that the k3 values start from 0.22. Does it ?

Hélène
POSTED BY: Hélène PASQUIER

Hi Jim, first of all thank you for your kind attention. Now I succeed. I've simplified my model a bit. K5 was not very useful ! I still have to implement this with my 10 kinetics...

Once again, many thanks! Hélène
POSTED BY: Hélène PASQUIER

Dear Jim The version I used is the one I found at the link you sent me, that is to say

Options[MultiNonlinearModelFit] = {AccuracyGoal -> Automatic, 
   ConfidenceLevel -> 19/20, EvaluationMonitor -> None, 
   Gradient -> Automatic, MaxIterations -> Automatic,
    Method -> Automatic, PrecisionGoal -> Automatic, 
   StepMonitor -> None, Tolerance -> Automatic, 
   VarianceEstimatorFunction -> Automatic,
    Weights -> Automatic, WorkingPrecision -> Automatic, 
   "DatasetIndexSymbol" -> \[FormalN]
   }; (* \[Equal] Join[Options[NonlinearModelFit], \
{"DatasetIndexSymbol" -> \[FormalN]}] *)

MultiNonlinearModelFit[datasets_, form_, fitParams_, 
   independents : Except[_List], opts : OptionsPattern[]] := 
      MultiNonlinearModelFit[datasets, form, 
   fitParams, {independents}, opts];

MultiNonlinearModelFit[datasets_, form : Except[{__Rule}, _List], 
   fitParams_, independents_, opts : OptionsPattern[]] := 
      MultiNonlinearModelFit[
   datasets, <|"Expressions" -> form, "Constraints" -> True|>, 
   fitParams, independents, opts];

MultiNonlinearModelFit[
       datasets : {__?(MatrixQ[#1, NumericQ] &)}, 
       KeyValuePattern[{"Expressions" -> expressions_List, 
     "Constraints" -> constraints_}], 
       fitParams_List, 
       independents_List,
       opts : OptionsPattern[]
   ] := Module[{
        fitfun, weights,
        numSets = Length[expressions], 
        augmentedData = Join @@ MapIndexed[
               Join[ConstantArray[N[#2], Length[#1]], #1, 2] &,
               datasets
           ], 
        indexSymbol = OptionValue["DatasetIndexSymbol"]
    },
       Condition[
            fitfun = With[{

       conditions = 
        Join @@ ({#1, expressions[[#1]]} & ) /@ Range[numSets]
               }, 
                  Switch @@ Prepend[conditions, Round[indexSymbol]]
              ]; 
            weights = Replace[
                  OptionValue[Weights],
                  {
                       (list_List)?(VectorQ[#1, NumericQ] & ) /; 
         Length[list] === numSets :> 

        Join @@ MapThread[ConstantArray, {list, Length /@ datasets}], 

       list : {__?(VectorQ[#1, NumericQ] & )} /; 
         Length /@ list === Length /@ datasets :>
                            Join @@ list, 

       "InverseLengthWeights" :> 
        Join @@ (ConstantArray[1./#1, #1] & ) /@ Length /@ datasets
                   }
              ]; 
            NonlinearModelFit[
                 augmentedData,

     If[TrueQ[constraints], fitfun, {fitfun, constraints}], 
                 fitParams,
                 Flatten[{indexSymbol, independents}],
                 Weights -> weights, 

     Sequence @@ FilterRules[{opts}, Options[NonlinearModelFit]]
             ]
            ,
            numSets === Length[datasets]
        ]
   ];

Thanks for the typo, I had forgotten to correct it in my message. Hélène
POSTED BY: Hélène PASQUIER
Jim Baldwin
Jim Baldwin, Retired
Posted 4 years ago

What version do you have that has a function called MultiNonlinearModelFit?

Found it: https://resources.wolframcloud.com/FunctionRepository/resources/MultiNonlinearModelFit.

And you do have a typo:

{agFunc, goxFunc, hFunc, fFunc} = {ag, gox, h, f} /. sol;

should be

{agFunc, goxFunc, hFunc, fFunc} = {ag, gox, h, ff} /. sol;
POSTED BY: Jim Baldwin
Reply to this discussion
Community posts can be styled and formatted using the Markdown syntax.
Tag limit exceeded
Note: Only the first five people you tag will receive an email notification; the other tagged names will appear as links to their profiles.
Reply Preview
Attachments
Remove
or Discard

Group Abstract Group Abstract