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/