MathGroup Archive 2001

[Date Index] [Thread Index] [Author Index]

Search the Archive

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/




  • Prev by Date: Re: Factorising operators??
  • Next by Date: Re: Q: Solving econometric models
  • Previous by thread: ListDensityPlot print problem + poll
  • Next by thread: Legend in ParametricPlot