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>
- NonlinearRegress and numerical functions...