MathGroup Archive 2006

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

Search the Archive

Re: FindFit and NormFunction

  • To: mathgroup at smc.vnet.net
  • Subject: [mg64217] Re: FindFit and NormFunction
  • From: Maxim <m.r at inbox.ru>
  • Date: Mon, 6 Feb 2006 02:49:19 -0500 (EST)
  • Sender: owner-wri-mathgroup at wolfram.com

-----Original Message-----
From: "David W. Cantrell" <DWCantrell at sigmaxi.org>
To: mathgroup at smc.vnet.net
Subject: [mg64217] Re: FindFit and NormFunction

> 
> Maxim <m.r at inbox.ru> wrote:
> > On Fri, 3 Feb 2006 00:11:19 +0000 (UTC), Paul Abbott
> > <paul at physics.uwa.edu.au> wrote:
> >
> > > I am having difficulty using the NormFunction option to FindFit. Let me
> > > give a concrete example. Abramowitz and Stegun Section 17.3.35 gives a
> > > nonlinear approximant to the complete elliptic integral. See
> > >
> > >  www.convertit.com/Go/ConvertIt/Reference/AMS55.ASP?Res=150&Page=592
> > >
> > > This approximant to EllipticE[1-m] can be implemented as
> > >
> > >  f[a_, b_, c_, d_][m_] = a m + b m^2 - (c m + d m^2) Log[m] + 1;
> > >
> > > After sampling EllipticE[1-m],
> > >
> > >  data = N[Table[{m, EllipticE[1-m]}, {m, 10^-8, 1 - 10^-8, 10^-3}]];
> > >
> > > using FindFit gives quite a decent approximant:
> > >
> > >  FindFit[data, f[a, b, c, d][m], {a,b,c,d}, m]
> > >
> > >  best[m_] =f[a, b, c, d][m] /. %
> > >
> > > The maximal absolute fractional error is ~6 x 10^-5 as seen from
> > >
> > >  Plot[10^5 (1 - best[m]/EllipticE[1 - m]), {m, 0, 1}, PlotRange -> All]
> > >
> > > However, the Abramowitz and Stegun approximant has error ~4 x 10^-5.
> > >
> > >  AS[m_] =f[0.4630151, 0.1077812, 0.2452727, 0.0412496][m];
> [snip]
> >
> > One option is to try a global optimizer instead of FindFit:
> >
> > In[4]:= (sol = NMinimize[
> >    Norm[data[[All, 2]] - f[a, b, c, d][data[[All, 1]]], Infinity],
> >    {a, b, c, d},
> >    Method -> {DifferentialEvolution, CrossProbability -> .05}]) // Timing
> >
> > Out[4]= {22.953*Second, {0.000035864209, {a -> 0.46224543, b ->
> > 0.10851758, c -> 0.24549621, d -> 0.041990523}}}
> 
> I've been unable to reproduce your result on my computer using version 5.1:
> 
> In[11]:= (sol = NMinimize[
>  Norm[data[[All,2]] - f[a, b, c, d][data[[All,1]]], Infinity],
>  {a, b, c, d},
>  Method -> {DifferentialEvolution, CrossProbability -> 0.05}]) // Timing
> 
> Out[11]= {14.422*Second, {0.000170617785, {a -> 0.47491173, b ->
> 0.09573489, c -> 0.24082805, d -> 0.034393029}}}
> 
> I also tried using several different values of RandomSeed, but was never
> able to get an error as low as your 0.000035864209. Why was I unable to
> reproduce your result? Could it be because I'm using version 5.1?
> 
> A triviality: Looking at the documentation for NMinimize, I had thought
> that double quotation marks would be _required_ around the Method option
> names, that is, Method -> {"DifferentialEvolution", "CrossProbability" ->
> 0.05}. But I see now that the quotation marks are not needed (although they
> still can be used). Is the documentation out of date?
> 
> > Plot[
> >    {EllipticE[1 - m] - AS[m],
> >     EllipticE[1 - m] - f[a, b, c, d][m] /. sol[[2]]},
> >    {m, 0, 1}, PlotRange -> All, PlotStyle -> {Red, Blue}]
> >
> > Actually this approximation has a smaller relative error as well as
> > smaller absolute error than AS[m].
> 
> Indeed it does. Explanation:
> 
> Formula 17.3.35 in A&S is from a book of computer approximations by
> Hastings. I haven't looked at that book, and so I don't know whether or not
> Hastings makes it clear that 1 + a + b must be Pi/2. But whether stated or
> not, that was apparently a restriction used by Hastings to insure that the
> approximation would give the exact value of EllipticE[1 - m] at the end
> point m = 1. And of course, the fact that the constant term in the
> approximation was 1 insured that it would also give the exact value at the
> other end point, m = 0. I assume that Hastings _wanted_ the approximation
> to be exact at both extremes of m. If my assumption is correct, then the
> values of the coefficients and the error given by Hastings are OK.
> 
> OTOH, if the only objective were to minimize worst |error|, then your
> result (which does not give the exact value at m = 1) is, as you indicated,
> better than Hastings. Furthermore, there is no reason that the constant
> term should be 1 if we do not require an exact value at m = 0. Dropping
> that restriction, we can then get a still better approximation,
> having |error| < 3*10^-5:
> 
>  f[a_, b_, c_, d_, e_][m_] = a m + b m^2 - (c m + d m^2) Log[m] + e
> 
> with {a -> 0.464167, b -> 0.106570, c -> 0.244708, d -> 0.040691,
> e -> 1.00002943}
> 
> One could similarly "improve" formula 17.3.36 in A&S. But to my taste at
> least, it's nicer to have approximations which give EllipticE[1 - m]
> exactly at the extremes of m.
> 
> Regards,
> David Cantrell
> 

You're correct, the NMinimize algorithms, in particular the local search methods (e.g., submethod InteriorPoint for RandomSearch), underwent some changes in version 5.2. My guess is that if you were to try

In[4]:= NMinimize[
  Norm[data[[All, 2]] - f[a, b, c, d][data[[All, 1]]], Infinity],
  {a, b, c, d},
  Method -> {DifferentialEvolution,
    CrossProbability -> .05, PostProcess -> False},
  MaxIterations -> 1000]

Out[4]= {0.000035461782, {a -> 0.46215121, b -> 0.10860872, c -> 0.24552379, d -> 0.042078762}}

(turning the local search off), the results would be the same in versions 5.1 and 5.2. Or we can optimize the function f[a, Pi/2 - 1 - a, c, d][m] if the value at m = 1 is fixed.

It is probably safer to use strings, as shown in the documentation ("DifferentialEvolution", etc.). For instance, I think in earlier versions of Mathematica Developer`SetSystemOptions accepted only strings, not symbols; now it accepts both. I just prefer to use symbols for uniformity and also because then I can use Ctrl-K (command completion) to speed up the typing.

Maxim Rytin
m.r at inbox.ru


  • Prev by Date: Re: notebook font encoding
  • Next by Date: thanks
  • Previous by thread: Re: FindFit and NormFunction
  • Next by thread: (Off-Topic): Re: Re: FindFit and NormFunction