# How do I find the best parameter to fit my data if the model is a interpolating function?

Hi I have a question regarding to find the best parameters for my model to fit my data. I have 3 ordinary equation, and I now just picked some parameters (k1 = 7.32*10^-5; k2 = 8.09*10^-9) for my ODE.

Defined all the constant values:

k1 = 7.32*10^-5;
k2 = 8.09*10^-9;
Da = 1.33*10 - 4;
Nd0 = 0.001;
HR0 = 0.2;
H0908 = 0.004;
x0 = 0.007;


And define the 3 ODE equations:

Equation1908 = {Nd0*
Caf908'[t] == -k1*(Nd0*HR0)/H0908*Caf908[t]*Cbf908[t]/Hf908[t] +
k2*(0.2 - 3*HR0*Cbf908[t]),
HR0*Cbf908'[t] == -3*k1*(Nd0*HR0)/H0908*Caf908[t]*
Cbf908[t]/Hf908[t] + 3*k2*(0.2 - 3*HR0*Cbf908[t]),
H0908*Hf908'[t] ==
3*k1*(Nd0*HR0)/H0908*Caf908[t]*Cbf908[t]/Hf908[t] -
3*k2*(0.2 - 3*HR0*Cbf908[t]),
Caf908 == 1, Cbf908 == 1, Hf908 == 1};


Using NDSolve:

BC1908 = NDSolve[
Equation1908, {Caf908[t], Cbf908[t], Hf908[t]}, {t, 0, 1000}]


At the end, I obtained 3 interpolating function, for Caf908[t], Cbf908[t], and Hf908[t]... I used the solution (Caf908[t]) to fit my data. X-axis (tfexp0908) versus y-axis (fexp0908):

fexp0908 = .001*{1.009790523, 0.898335138, 0.878948419,
0.830114856, 0.767123385, 0.732170062, 0.672106602,
0.637589428, 0.59141947, 0.523944512, 0.554584169,
0.451444203, 0.396545111, 0.444125908, 0.352355452,
0.272913948, 0.33877861, 0.287900412, 0.276936425};
tfexp0908 = {0, 15, 30, 45, 60, 75, 90, 105, 120, 135, 150, 180, 210,
245, 270, 300, 330, 360, 390};


And I obtained a good fit. But I would like to know how I can use Mathematica to obtain the best parameter, k1 & k2, for my model (Caf908[t]) and get the best fit to my data. Thank you!

• its doable in principle but right now your set of equations does not seem to integrate properly. You forgot to give k1 a value Nov 4, 2012 at 20:23
• Hi Chris, Thanks for your reply. I'm sorry that I forgot to include a value for k1. K1& k2 values are:k1 = 7.32*10^-5; k2 = 8.09*10^-9; Nov 4, 2012 at 21:46

Let's clear the k1 and k2 variables and put in kk1 and kk2 the corresponding values

Clear[k1, k2];
kk1 = 7.32*10^-5;
kk2 = 8.09*10^-9;


Now define a set of equations to solve depending on k1,k2

Equation[k1_?NumberQ,k2__?NumberQ] = {Nd0*
Caf908'[t] == -k1*(Nd0*HR0)/H0908*Caf908[t]*Cbf908[t]/Hf908[t] +
k2*(0.2 - 3*HR0*Cbf908[t]),
HR0*Cbf908'[t] == -3*k1*(Nd0*HR0)/H0908*Caf908[t]*
Cbf908[t]/Hf908[t] + 3*k2*(0.2 - 3*HR0*Cbf908[t]),
H0908*Hf908'[t] ==
3*k1*(Nd0*HR0)/H0908*Caf908[t]*Cbf908[t]/Hf908[t] -
3*k2*(0.2 - 3*HR0*Cbf908[t]), Caf908 == 1, Cbf908 == 1,
Hf908 == 1};


Now let us define a $\chi^2$ function

Clear[chi2];
chi2[k1_?NumberQ, k2_?NumberQ] :=
(sol =
Caf908 /. NDSolve[
Equation[k1, k2], {Caf908, Cbf908, Hf908}, {t, 0, 1000}] //
First;
(sol /@ tfexp0908) - fexp0908 // #.# &
)


and let's sample this function around what is supposedly the correct value

Table[chi2[kk1*i, kk2*j], {i, 0.1, 2.1, 0.1}, {j, 0.1, 5.1,
0.2}] // ListContourPlot It seems to suggest i) the guess is not so good ii) somehow the $\chi^2$ does not depend on k2

Finding the right range of values for kk1 and kk2 and sorting out why the $\chi^2$ does not depend on k1 is beyond a Mathematica question?

Note that one could use in principle NMinimize to find the solution more efficiently than via sampling the likelihood function, but before these points are addressed there's no point.

EDIT

I believe you should remove the 0.01* in front of fexp0908

Clear[chi2];
chi2[k1_?NumberQ, k2_?NumberQ,
debug_: False] := (sol =
Caf908 /.
NDSolve[Equation[k1, k2], {Caf908, Cbf908, Hf908}, {t, 0,
1000}] // First;
If[debug, {(sol /@ tfexp0908), fexp0908} // ListLinePlot // Print];
(sol /@ tfexp0908) - fexp0908 // #.# &)


because without it I get

 Table[{i*kk1,chi2[kk1*i, kk2*j]},{i, 0.9, 2.3, 0.1},
{j, 0.1, 5.1, 0.2}]//Transpose//ListLinePlot and Indeed the fit seems okish:

 chi2[kk1*1.3, kk2, True] Now the minimization

NMinimize[{chi2[k1, k2], k1 > 0, k2 > 0}, {k1, k2}]

(* {0.0339566,{k1->0.000110757,k2->0.} *)


yields a poorer fit but hasn't converged.

It seems your set of data poorly constrains k2

• Thanks for your help, Chris. Do you know if there is anyway I can find k1&k2 if the values are incorrect? I just started learning mathematica, and know there are findfit and fit. But I must have a function in order to use these two method, and I only have interpolating function from NDSolve. Could you show me how to use NMinimize? Thanks! Nov 4, 2012 at 21:50
• Hi Chris, the problem is for solving chemical reaction, and k1&k2 are the reaction constants of forward and reverse reactions. I'm still having some problems though. When I use k1&k2 (k1->15.3735,k2->2.45184*10^-9) obtained from NMinimize to fit my data, it gave me a very bad fit, while the previous values (k1 = 7.32*10^-5; k2 = 8.09*10^-9) fits better still. Thanks. Nov 5, 2012 at 0:04
• Thank you so much, Chris! I finally solved the problem after fixing some values and my equations!! I'll let you know if I have more questions! Nov 5, 2012 at 20:13
• @chris What is the //#.#&` at the end of the chi2 function? Mar 21, 2013 at 18:45
• NVM. I just realized it's the square of the residuals between what's predicted by the interpolation and the actual data. Mar 21, 2013 at 19:03