1
$\begingroup$

How does one use a fitting routine such as NonlinearModelFit with a ParametricFunction? Example below:

First let's generate some data:

sample[t_] = (0.002 + 101 t - 461000 t^2 + 2.218 10^9 t^3 - 
             3.64 10^12 t^4 + 3.17 10^15 t^5 ) Exp[-8653 t];
data = Table[{t, sample[t] + 
       RandomVariate[NormalDistribution[0, 0.00001]]}, 
       {t, 0, 0.002, 0.000004}];
ListPlot[data]

enter image description here

Now lets define a model:

rateeqs = {a'[t] == -k1a a[t] - k12 a[t] c70gs + k21 b[t] c60gs, 
           b'[t] == -k1b b[t] - k21 b[t] c60gs + k12 a[t] c70gs, 
           a[0] == a0, b[0] == b0};
c60gs = c70gs = 5;
maxTime = 0.0025;
e60 = 19060;
e70 = 948;
fitFunc[t_] = e60 a[t] + e70 b[t];
params = {k1a, k1b, k12, k21, a0, b0};
initGuesses = {8000, 100, 4500, 2000, 5. 10^-8, 8 10^-7};

Now we generate a parametric solution to the model:

solution[t_] = ParametricNDSolveValue[rateeqs, fitFunc[t], {t, 0, maxTime}, 
               params, WorkingPrecision->30]

And we verify that the solution works:

Show[ListPlot[data, PlotStyle -> Red, PlotRange -> Full], 
     Plot[(solution[t] @@ initGuesses), {t, 0, maxTime}]]

enter image description here

It would seem that I should be able to fit this model to the sample data starting from the initGuesses values for the parameters using something like this:

modelFit = NonlinearModelFit[data, solution[t] @@ params, 
           {params, initGuesses}\[Transpose], t];

But I get errors:

enter image description here

I'm sure it's a syntax issue, but I have tried so many different possibilities that the screen is swimming in front of my eyes. Any help is appreciated!

Improve this question
$\endgroup$
2
  • 1
    $\begingroup$ You have typo in your post. You define $e60$ and $e70$, but they occur as greek $\epsilon60$ and $\epsilon70$ in your fitFunction[]. $\endgroup$
    – Tim Laska
    May 11, 2019 at 11:49
  • $\begingroup$ Thank you. As I was copying it in I tried to simplify from greek to e, but obviously missed a couple spots. I will correct. $\endgroup$
    – Kevin Ausman
    May 11, 2019 at 18:26

2 Answers 2

Reset to default
4
$\begingroup$

For some reason I can't get it to work with the new fix for solution, but it works with the workaround from the answer to your previous question.

solution = ParametricNDSolveValue[Join[rateeqs, {y'[t] == fitFunc'[t], y[0] == fitFunc[0]}], y, {t, 0, maxTime}, params];

fit = NonlinearModelFit[data, (solution @@ params)[t], params, t];

Show[
 ListPlot[data, PlotStyle -> Red, PlotRange -> Full],
 Plot[fit[t], {t, 0, maxTime}, PlotStyle -> Thick, PlotRange -> Full]]

enter image description here

fit["ParameterTable"]

enter image description here

Edited to include parameter guesses

Sometimes re-running this yields different (worse) results (presumably due to different samples for data), and a warning about getting stuck in a local minimum.

To provide starting values:

newGuesses = {1000, 30000, 400, -100, 0, 0};
fit2 = NonlinearModelFit[data, (solution @@ params)[t], 
  Transpose[{params, newGuesses}], t]
Improve this answer
$\endgroup$
2
  • $\begingroup$ Hmmm. I was thinking the issue was the substituting in of [t], and one advantage of your approach is that the output expression doesn’t have an explicit [t]. Thank you for the suggestions; I will give this a try. $\endgroup$
    – Kevin Ausman
    May 11, 2019 at 3:15
  • $\begingroup$ And yes, this works. Thank you very much! $\endgroup$
    – Kevin Ausman
    May 11, 2019 at 18:27
2
$\begingroup$

I recently answered another chemical kinetics fitting to experimental data question Here. It uses Manipulate to populate the initial parameter guesses for the fit. I also try to mimic some functionality from the nonlinear curve fitting platform of the statistical package JMP using a string template. JMP provides a Manipulate like slider model to adjust parameters post fit. Real world data will not be as well behaved as synthetic data.

Setting Up Slider Model For Initial Guesses

I modified your approach as to be more consistent with my previous posted approach. Now, we will create a slider model with global variable to track the parameters with Manipulate. I also plotted the underlying basis functions to get a better idea which ones are driving over model features.

(* Synthesized Data *)
sample[t_] = (0.002 + 101 t - 461000 t^2 + 2.218 10^9 t^3 - 
     3.64 10^12 t^4 + 3.17 10^15 t^5) Exp[-8653 t];
data = Table[{t, 
    sample[t] + RandomVariate[NormalDistribution[0, 0.00001]]}, {t, 0,
     0.002, 0.000004}];
(*Parameters*)
c60gs = c70gs = 5;
maxTime = 0.0025;
e60 = 19060;
e70 = 948;
(*Equations*)
eqnA = a'[t] == -k1a a[t] - k12 a[t] c70gs + k21 b[t] c60gs;
eqnB = b'[t] == -k1b b[t] - k21 b[t] c60gs + k12 a[t] c70gs;
ics = {a[0] == a0, b[0] == b0};
eqns = {eqnA, eqnB}~Join~ics;
(*Variables*)
vbles = {a, b};
(*Parameters*)
params = {k1a, k1b, k12, k21, a0, b0};
pfun = ParametricNDSolveValue[eqns, vbles, {t, 0, maxTime}, params, 
  MaxStepFraction -> 0.001]
(*Create an appropriate model function to fit*)
model[k1a_, k1b_, k12_, k21_, a0_, b0_][
   t_] := {e60 #[[1]], e70 #[[2]], e60 #[[1]] + e70 #[[2]]} &@
    Through[pfun[k1a, k1b, k12, k21, a0, b0][t], List] /; 
   And @@ NumericQ /@ {k1a, k1b, k12, k21, a0, b0};
(*Slider Model to Set Initial Parameter Guesses*)
Manipulate[global = {k1a, k1b, k12, k21, a0, b0};
 Show[ListPlot[data, PlotStyle -> Red, PlotRange -> Full, 
   PlotLegends -> {"data"}], 
  Plot[Evaluate@model[k1a, k1b, k12, k21, a0, b0][t], {t, 0, maxTime},
    PlotRange -> Full, 
   PlotLegends -> {"a[t]", "b[t]", "model[t]"}]], {{k1a, 8000}, 4000, 
  15000, Appearance -> "Labeled"}, {{k1b, 100}, 0, 800, 
  Appearance -> "Labeled"}, {{k12, 4500}, 2000, 6500, 
  Appearance -> "Labeled"}, {{k21, 2000}, 1000, 4000, 
  Appearance -> "Labeled"}, {{a0, 5. 10^-8}, 0, 10 5. 10^-8, 
  Appearance -> "Labeled"}, {{b0, 8 10^-7}, 0, 10 8 10^-7, 
  Appearance -> "Labeled"}]
(*Display global variable*)
Dynamic@global

Initial Manipulate

Setting Up The String Template

The Linux world has the concept of a "Here Document", which makes it quite easy to take an existing piece of code and make it a template. The Mathematica equivalent is the StringTemplate[] or TemplateApply[] functions. We will use it to replace the initial values of the Manipulate function to represent what was returned by the fit. We can delimit template variables with back ticks "`". I chose this approach because I was using a similar looking symbol for many different uses (replacement rule, manipulate variables, etc). There could be a better way, but this worked. Using the manipulate from above (remembering to escape the quotes), I came up with the following string template:

strtemp = 
  "Manipulate[Show[ListPlot[data,PlotStyle->Red,PlotRange->Full,\
PlotLegends->{\"data\"}],Plot[Evaluate@model[k1a,k1b,k12,k21,a0,b0][t]\
,{t,0,maxTime},PlotRange->Full,PlotLegends->{\"a[t]\",\"b[t]\",\"\
model[t]\"}]],{{k1a,`k1a`},4000,15000,Appearance->\"Labeled\"},{{k1b,`\
k1b`},0,800,Appearance->\"Labeled\"},{{k12,`k12`},2000,6500,\
Appearance->\"Labeled\"},{{k21,`k21`},1000,4000,Appearance->\"Labeled\
\"},{{a0,`a0`},0,10 5. 10^-8,Appearance->\"Labeled\"},{{b0,`b0`},0,10 \
8 10^-7,Appearance->\"Labeled\"}]";

Fitting the Parametric Equation

We populate the our initial parameter guesses using the "global" and Dynamic from our previous Manipulate to drive a fit and then create a new Manipulate initialized with the fitted parameters using the string template as shown below. 

initguess = 
  MapThread[
   List, {{k1a, k1b, k12, k21, a0, b0}, First@Dynamic@global}];
fit = FindFit[data, model[k1a, k1b, k12, k21, a0, b0][t][[3]], 
  initguess, t, Method -> NMinimize]
ToExpression@
 TemplateApply[strtemp, 
  AssociationThread @@ ({ToString@# & /@ Keys@fit, Values@fit})]

Fitted Image

Now, you have the ability to make small adjustments to the model to account for things such as experimental artifacts that are bound to occur in real life.

Improve this answer
$\endgroup$
2
  • $\begingroup$ Thank you for this. I already have developed almost all of this type of functionality on my own, plus a number of other valuable features (graphical parameter sensitivity analysis, automatic generation of differential rate laws from arbitrary chemical mechanisms where the mechanisms are entered in standard chemistry notation, and several others), and am just asking about the few bits that were missing. So don't take my lack of using your code as dismissal; your work looks great. I just am not looking to reinvent the wheel I've already invented. $\endgroup$
    – Kevin Ausman
    May 11, 2019 at 18:35
  • $\begingroup$ @KevinAusman No problem. I just wanted to show a workflow to the community how Mathematica could be used to generate better initial guesses efficiently and to create a useful JMP like functionality to make minor adjustments to fitted parameters. Also, I was uneasy with the current accepted answer generating unphysical parameter estimates for parameters that should be intrinsically positive. I thought this workflow would be more robust by seeding the fit with better initial parameter estimates. $\endgroup$
    – Tim Laska
    May 11, 2019 at 20:01

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.