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:

> >
> > 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
> >
> >
> > _________________________________________________
> > 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}
>
>
> 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

--

```

• Prev by Date: Re: Why Doesn't N[Pi,i] Give i Digits For Small i, Mathematica
• Next by Date: Re: Help Needed: How to use a Notebook from an external Program via Mathlink?
• Previous by thread: Re: Re: NonlinearRegress and numerical functions...
• Next by thread: RE: Re: NonlinearRegress and numerical functions...