4

This is my first post, and I will preface my question by stating that I am relatively new to Mathematica.

In any case, I would like to fit the solution of chemical kinetics equations to data and get the resulting "MeanPredictionBands". Using the answer to a previous question as a starting point, I am able to fit a very simple set of chemical kinetics equations to two data sets self-consistently (see example code below).

tmax = 100.; (*region of time in which we are interested to have a \
solution*)

(*Define a simple chemical kinetics coupled differential equation*)
sol = ParametricNDSolveValue[
  {D[s1[t], t] == k1 s1[t] s2[t],
   D[s2[t], t] == -k1 s1[t] s2[t],
   s1[0] == 0.1,
   s2[0] == 10.},

  {s1, s2},

  {t, 0, tmax},

  {k1}
  ];

(*temp variables that help us avoid solving the coupled differential \
equations over and over if "model" is called without changing k1temp*)


k1temp = 0.1;
soltemp = sol[k1temp];

(*the function definition that will return the solution of the \
coupled differential equations for a given k1, species and time*)
model[k1_?NumericQ][species_?NumericQ, t_?NumericQ] :=
  Module[{returnval, intspecies},
   If[k1 != k1temp,
    soltemp = sol[k1];
    k1temp = k1;
    ];
   intspecies = Round[species];
   soltemp[[intspecies]][t]
   ];

(*Make some fake data using model*)
data = {};
Do[
  Do[
   noise = RandomReal[{-0.5, 0.5}];
   AppendTo[data, {i, t, model[0.2][i, t] + noise}];
   , {t, 0., tmax/10.}]
  , {i, 1, 2}];

(*fit the fake data*)
fit = NonlinearModelFit[data, model[k1f][i, t], {{k1f, 0.3}}, {i, t}];

fit["ParameterTable"]

Using the FittedModel from NonlinearModelFit, I expected it to be straight forward to get the "MeanPredictionBands" via the following function definition.

meanpredictionbandsfunc[i_?NumericQ, t_?NumericQ] = 
 fit["MeanPredictionBands"]

This function definition, however, is giving me a number of warnings/errors and when evaluated it does not give me the "MeanPredictionBands" as I would expect. For example, the evaluation

meanpredictionbandsfunc[1, 10]

is giving me the following errors

StringForm::sfr: Item 2 requested in "`1` cannot be interpreted. The operator `2` requires a subscript with a variable specification." out of range; 1 items available.

Internal`LocalizedBlock::novar: 1 cannot be interpreted. The operator `2` requires a subscript with a variable specification.

StringForm::sfr: Item 2 requested in "`1` cannot be interpreted. The operator `2` requires a subscript with a variable specification." out of range; 1 items available.

Internal`LocalizedBlock::novar: 1 cannot be interpreted. The operator `2` requires a subscript with a variable specification.

StringForm::sfr: Item 2 requested in "`1` cannot be interpreted. The operator `2` requires a subscript with a variable specification." out of range; 1 items available.

Internal`LocalizedBlock::novar: 1 cannot be interpreted. The operator `2` requires a subscript with a variable specification.

How should I define a function to provide me with the "MeanPredictionBands"?

I would really appreciate your help/suggestions!

Community
  • 1
Rischaard
  • 41
  • 3

1 Answers1

3

In this case because your model is not simple (in that it contains If statements and interpolating functions) both MeanPredictionBands and SinglePredictionBands don't work. (I think that use of functions like Abs in the model also makes MeanPredictionBands not work - but I could be wrong about that.)

I've removed the If statements from your model but I think it works identically as before:

model[k1_?NumericQ, species_?NumericQ, t_?NumericQ] := sol[k1][[species]][t];
fit = NonlinearModelFit[data, model[k1f, i, t], {{k1f, 0.3}}, {i, t}];

There are two things one can do. First, the values of the MeanPredictionBands can be obtained for the individual observations using MeanPredictionConfidenceIntervals:

fit["MeanPredictionConfidenceIntervals"]
(* {{0.1, 0.1}, {0.627004, 0.705878}, {3.07832, 3.64677}, {7.53791, 8.1988},
 {9.61518, 9.80514}, {10.0249, 10.0609}, {10.0888, 10.095}, {10.0983, 10.0994},
 {10.0998, 10.0999}, {10.1, 10.1}, {10.1, 10.1}, {10., 10.}, {9.39412, 9.473},
 {6.45323, 7.02168}, {1.9012, 2.56209}, {0.294855, 0.484816}, {0.0390947,  0.0750506},
 {0.00503452, 0.0112016}, {0.000639995, 0.00165984}, {0.0000802901,  0.000245284},
 {9.9049*10^-6,  0.0000361818}, {1.19425*10^-6, 5.32897*10^-6}} *)

The other thing to do is to construct the confidence bands directly using the CovarianceMatrix and the Delta Method. For your example you only have one parameter so you can find the lower and upper bands just by plugging in the lower and upper confidence limits for k1f:

k1mle = k1f /. fit["BestFitParameters"];
{k1lower, k1upper} = fit["ParameterConfidenceIntervals"][[1]]

Show[ListPlot[data[[All, {2, 3}]]],
 Plot[{model[k1mle, 1, t], model[k1lower, 1, t], model[k1upper, 1, t]}, 
   {t, 0, Max[data[[All, 2]]]}, PlotStyle -> {Black, Green, Green}],
 Plot[{model[k1mle, 2, t], model[k1lower, 2, t], model[k1upper, 2, t]},
   {t, 0, Max[data[[All, 2]]]}, PlotStyle -> {Black, Green, Green}]]

Confidence bands

There are two bands because you have different curves depending on species.

JimB
  • 41,653
  • 3
  • 48
  • 106
  • The “Delta Method” is straight forward for this 1-parameter, simple example I provided above, but what I really need is a solution when there are 5-10 parameters. Though, I admit, that the “Delta Method” can probably be applied in this case as well, though it will certainly not be this straight forward. – Rischaard Dec 05 '16 at 17:42
  • I don’t think the “If” statements need to be removed if you want to use fit[“MeanPredictionConfidenceIntervals”]. Removing the “If” statements from the definition of "model", as you propose, results in the same value, but can be very slow when used in "Plot[]", especially when it takes a long time (i.e. minutes or more) to solve "sol" since “sol” is evaluated for each point in time and really only needs to be evaluated once for a given set of parameters. – Rischaard Dec 05 '16 at 17:43