Using NonlinearRegress with Dif. Equations
- To: mathgroup at smc.vnet.net
- Subject: [mg30979] Using NonlinearRegress with Dif. Equations
- From: guillerm at usal.es
- Date: Sat, 29 Sep 2001 04:18:58 -0400 (EDT)
- Sender: owner-wri-mathgroup at wolfram.com
(*A excellent method for fitting experimental data and differential equations
(using nonlinearRegress) was shown to mathgroup by Carl Woll. This method need
that the user teach to Mathematica the derivatives for each coeffs fo be
fitted. It can be complex when the number of variables and coefficients is
big. I have found that the problem can be made esier using the numeric
derivatives applying the package NumericalMath`NLimit`. I will show with an
example: The problem is to fit the coeffs. a12 and a23 in the bellow SODE using
experimental data given in list dataIodine *)
In[1]:= model[(t1_)?NumberQ, (a12_)?NumberQ, (a23_)?NumberQ] :=
x2[t] /. First[NDSolve[{Derivative[1][x1][t] == 34.4/E^(200.*t) + 1.09/E^
(110.*t) +
0.808/E^(102.*t) + 6.414/E^(100.*t) + 5.458/E^(24.*t) + (-1.9404 - a12)*x1[t] +
0.0462*x3[t], Derivative[1][x2][t] == a12*x1[t] - a23*x2[t],
Derivative[1][x3][t] == a23*x2[t] - 0.05775*x3[t], x1[0] == 0, x2[0] == 0,
x3[0] == 0}, {x1[t], x2[t], x3[t]}, {t, 0, 500}, {AccuracyGoal -> 100,
MaxSteps -> 3000}]] /. t -> t1
In[2]:=
dataIodine = {{0, 0.002}, {50, 0.106}, {100, 0.077}, {150, 0.056}, {200,
0.041}, {250, 0.032}, {300, 0.023}, {350, 0.018}, {400, 0.011}, {450, 0.011},
{500, 0.007}};
(*Because I want to apply the gradient method I teach the derivative to
Mathematica.Usually it is computer the derivatives using the de differential
equations but I prefer this way that is easier*)
In[3]:=
Needs["NumericalMath`NLimit`"]
In[4]:=
Derivative[1,0,0][model][n_,a1_,a3_]:=ND[model[a1,a2,a3],a1,n]
(*In this example we don?t need this derivative, but it is included to show the
pattern *)
In[5]:=Derivative[0,1,0][model][a1_,n_,a3_]:=ND[model[a1,a2,a3],a2,n]
In[6]:=Derivative[0,0,1][model][a1_,a2_,n_]:=ND[model[a1,a2,a3],a3,n]
You can realize that this pattern could be apply to others model with two
parameters for fitting. It will taken a long time but the solution look like
fine (I have apply this method for other more complex model with 5 coefficients
to be fitted and it works)
In[7]:=
Needs["Statistics`NonlinearFit`"];
NonlinearRegress[dataIodine,model[t,a12,a23],{t},{{a12,{0.8,1}},{a23,
{0.007,0.008}}}, Method\[Rule] Gradient,RegressionReport\[Rule]
{BestFit,EstimatedVariance,ParameterCITable, ANOVATable,FitResiduals}]//Timing
(*Out[8]:=
{1299.54*Second, {BestFit -> model[t, 0.8657361750616687, 0.00842616513059947],
EstimatedVariance -> 1.6307890484089922*^-6, ParameterCITable ->
TableForm[{{"", "Estimate", "Asymptotic SE", "CI"}, {a12,
0.8657361750616687,
0.032891536986871514, {0.7913303490713752, 0.9401420010519621}},
{a23, 0.00842616513059947, 0.00045863543397038983, {0.007388659698530291,
0.00946367056266865}}}, TableDepth -> 2, TableHeadings ->
{{a12, a23}, {"Estimate", "Asymptotic SE", "CI"}}],
ANOVATable -> TableForm[{{2, 0.024139322898564318, 0.012069661449282159},
{9, 0.000014677101435680929, 1.6307890484089922*^-6}, {11, 0.024154},
{10, 0.010748909090909092}},.....*)
(*I am working in expanding this method building a function that using the
length of coefficients list to be fitted for generate the derivative.
ie: For three parameter for fitting the derivative for apply will be (Let
observer that will need the derivative for 4 parameter to be account the time t
Derivative[1,0,0,0][headmodel][n_,a2_,a3_,a4_]:=
ND[headmodel[a1,a2,a3,a4],a1,n];
Derivative[0,1,0,0][headmodel][a1_,n_,a3_,a4_]:=
ND[headmodel[a1,a2,a3,a4],a2,n];
Derivative[0,0,1,0][headmodel][a1_,a2_,n_,a4_]:=
ND[headmodel[a1,a2,a3,a4],a3,n];
Derivative[0,0,0,1][headmodel][a1_,a2_,a3_,n_]:=
ND[headmodel[a1,a2,a3,a4],a4,n];
and so on *)
(*I am thinking in build a function more or less like this. It wil be better if
the derivative are genereted as funtion of parameters to be fitted length, in
fact I have find solution but are not eleghant *)
In[8]:=
nonlinearRegressNumeric[data_, model_, variables_, parameters__] :=
Module[{n, a1, a2, a3, a4, headmodel}, headmodel = Head[model];
Derivative[1, 0, 0][headmodel][n_, a2_, a3_] := ND[headmodel[a1, a2, a3], a1,
n];
Derivative[0, 1, 0][headmodel][a1_, n_, a3_] := ND[headmodel[a1, a2, a3], a2,
n];
Derivative[0, 0, 1][headmodel][a1_, a2_, n_] := ND[headmodel[a1, a2, a3], a3,
n];
Derivative[1, 0, 0, 0][headmodel][n_, a2_, a3_, a4_] := ND[headmodel[a1, a2,
a3, a4], a1, n];
Derivative[0, 1, 0, 0][headmodel][a1_, n_, a3_, a4_] :=
ND[headmodel[a1, a2, a3, a4], a2, n]; Derivative[0, 0, 1, 0][headmodel][a1_,
a2_, n_,
a4_] := ND[headmodel[a1, a2, a3, a4], a3, n];
Derivative[0, 0, 0, 1][headmodel][a1_, a2_, a3_, n_] := ND[headmodel[a1, a2,
a3, a4], a4, n];
NonlinearRegress[data, model, variables, parameters, Method -> Gradient]]
In[9]:=
nonlinearRegressNumeric[dataIodine, model[t, a12, a23], {t}, {a12, {0.8, 1}},
{a23, {0.007, 0.008}}}, Method -> Gradient]
I will appreciate any help. If anybody has problems with this E-mail I can send
the Mathematica notebook
Guillermo
---------------------------------------------
This message was sent using Endymion MailMan.
http://www.endymion.com/products/mailman/