# Fitting experimental data with a parametric function

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]


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}]]


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:

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!

• You have typo in your post. You define $e60$ and $e70$, but they occur as greek $\epsilon60$ and $\epsilon70$ in your fitFunction[]. – Tim Laska May 11 at 11:49
• 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. – Kevin Ausman May 11 at 18:26

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]]



fit["ParameterTable"]


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]

• 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. – Kevin Ausman May 11 at 3:15
• And yes, this works. Thank you very much! – Kevin Ausman May 11 at 18:27

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


# 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 =
`