# 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[0] == 1, Cbf908[0] == 1, Hf908[0] == 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 – chris Nov 4 '12 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; – DumbleKo Nov 4 '12 at 21:46
• Welcome to Mathematica.SE! I suggest the following: 1) As you receive help, try to give it too, by answering questions in your area of expertise. 2)Read the FAQs! 3) When you see good Q&A, vote them up by clicking the gray triangles, because the credibility of the system is based on the reputation gained by users sharing their knowledge. ALSO, remember to accept the answer, if any, that solves your problem, by clicking the checkmark sign – chris Nov 5 '12 at 20:15

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[0] == 1, Cbf908[0] == 1,
Hf908[0] == 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! – DumbleKo Nov 4 '12 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. – DumbleKo Nov 5 '12 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! – DumbleKo Nov 5 '12 at 20:13
• @chris What is the //#.#& at the end of the chi2 function? – Christopher Bowman Mar 21 '13 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. – Christopher Bowman Mar 21 '13 at 19:03