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
- Follow-Ups:
- (Off-Topic): Re: Re: FindFit and NormFunction
- From: gardyloo <gardyloo@mail.wsu.edu>
- (Off-Topic): Re: Re: FindFit and NormFunction