Re: NonlinearRegress and numerical functions...

*To*: mathgroup at smc.vnet.net*Subject*: [mg20254] Re: [mg20161] NonlinearRegress and numerical functions...*From*: Mike Trefry <mike.trefry at per.clw.csiro.au>*Date*: Fri, 8 Oct 1999 18:30:25 -0400*Organization*: CSIRO Land and Water*References*: <199910040107.VAA17261@smc.vnet.net.> <7tcbbl$m6r@smc.vnet.net>*Sender*: owner-wri-mathgroup at wolfram.com

Hi Daniel, I think the essential point in this discussion is not that Mathematica won't find a minimum vector of parameters in a pde system, it is that Mathematica won't do a NonlinearRegress on a vector of parameters in a pde system. That is, we have no way of getting confidence intervals on the parameter values returned by things like FindMinimum, NonlinearFit, etc. This is a critical failure. Mike Daniel Lichtblau wrote: > Lars Ragnarsson wrote: > > > > Hi > > > > I've got a big problem making Mathematica realize that my function is > > numerical. The following problem illustrates it in a rather simplified > > way: > > > > g[t1_] := y[t] /. DSolve[{y'[t] == -y[t] 3, y[0] == 6}, y[t], t][[1]] /. > > > > t -> t1 > > > > data = {0, 10, 20, 100}; > > > > fit = Transpose[{data, g[data]}]; > > > > f[t1_, a_,b_] := (y[t] /.NDSolve[{y'[t] == -y[t] a, y[0] == b}, y[t], > > {t, 0, 100}][[1]]) /.t -> t1 > > > > NonlinearRegress[data, f[t, a, 6], {t}, {a, {2.9, 3.2}, 2, 4}] > > > > It should be rather straightforward to solve this, but NonlinearRegress > > > > seem to evaluate the function 'f' first and the results: > > > > NDSolve::"ndnum": "Encountered non-numerical value for a derivative at \ > > > > \!\(t\) == \!\(2.680466916142905`*^-274\)." > > > > etc.... > > > > I've tried to use Unevaluated and Hold but nothing works! > > > > Any help would be greatly appreciated! > > > > Regards > > > > Larske Ragnarsson > > > > _________________________________________________ > > Larske Ragnarsson > > Phone: +46 (0)31 7721867 > > Fax: +46 (0)31 7723622 > > Chalmers University of Technology > > Solid State Electronics Laboratory > > Department of Microelectronics ED > > S-412 96 Guteborg, SWEDEN > > _________________________________________________ > > I'll change the notation just a bit. To start we define our function, > create a bunch of data points, and show, with NonlinearFit, that we can > recover the function. > > g[a_, b_, t_] := (y[t1] /. > First[DSolve[{y'[t1]==-a*y[t1], y[0]==b}, y[t1], t1]]) /. t1->t > pnts = {0, 10, 20, 100}; > data = {#, g[3,6,t] /. t->#}& /@ pnts > <<Statistics` > > In[5]:= NonlinearFit[data, g[a,6,t], t, {a,1}] > 6 > Out[5]= ----- > 3. t > E > > Now suppose instead that we have a function defined by (numerical) > solutions to a differential equation. I assume the example given was for > purposes of comparison to a known solution though of course in general > we do not have this recourse. Let me point out some issues, at least one > of which is substantial. > > First, your numeric differential equation solution is not very accurate; > to obtain a result similar to the exact solution for x as large as 100 > you need to fiddle with the NDSolve options (maybe this numeric error > was intentional?). > > Clear[f] > f[a_?NumberQ, b_?NumberQ, t_?NumberQ] := > (y[t1] /. First[NDSolve[{y'[t1]==-a*y[t1], y[0]==b}, y[t1], > {t1,0,100}, AccuracyGoal->150, MaxSteps->2000]]) /. t1->t > > Let us check the relative errors at our data points. > > In[9]:= Map[(f[3,6,#]-g[3,6,#])/g[3,6,#]&, pnts] > Out[9]= {0., -0.0000162022, -0.0000385361, -0.000219838} > > Not bad. > > Now for a more serious issue. You want to do a nonlinear function fit to > a model that is, in effect, defined by a program. While this is a > reasonable thing to want, I am not certain it can be done directly in > Mathematica (I confess to no great expertise with NonlinearFit/Regress, > so I may be mistaken). Anyway, you can instead formulate this as a > minimization problem and use FindMinimum to solve for the parameter > value. We do this explicitly below. Again, some option fiddling is > helpful in order to get a good result. > > In[15]:= objfunc = Apply[Plus, Map[(f[a, 6, #[[1]]] - #[[2]])^2 &, > data]]; > > In[16]:= FindMinimum[objfunc, {a,2,2.1}, AccuracyGoal->50, > MaxIterations->200, WorkingPrecision->100] > -53 > Out[16]= {3.79385 10 , {a -> 2.99999837977173152356769137355}} > > Given the disparities between f and g at the data points, this is likely > about as accurate a solution as can be obtained; to do better one would > have to get more accurate results from the model function f and perhaps > also use/request higher accuracy/precision in FindMinimum. > > Daniel Lichtblau > Wolfram Research --

**References**:**NonlinearRegress and numerical functions...***From:*Lars Ragnarsson <loke@ic.chalmers.se>