Re: FindFit and NormFunction
- To: mathgroup at smc.vnet.net
- Subject: [mg64171] Re: FindFit and NormFunction
- From: Maxim <m.r at inbox.ru>
- Date: Fri, 3 Feb 2006 01:04:00 -0500 (EST)
- References: <dru737$a0v$1@smc.vnet.net>
- Sender: owner-wri-mathgroup at wolfram.com
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];
>
> Plot[10^5 (1 - AS[m]/EllipticE[1 - m]), {m, 0, 1}, PlotRange -> All]
>
> Essentially, one needs to fit with respect to the L-Infinity-Norm, see
>
> http://mathworld.wolfram.com/L-Infinity-Norm.html
>
> as the NormFunction in FindFit. However, when I try
>
> FindFit[data, f[a, b, c, d][m], {a,b,c,d}, m,
> NormFunction :> (Norm[#, Infinity]&)]
>
> I get a FindFit::"lmnl" message (which is reasonable and informative)
> and also a FindFit::"lstol" message:
>
> "The line search decreased the step size to within tolerance specified
> by AccuracyGoal and PrecisionGoal but was unable to find a sufficient
> decrease in the norm of the residual. You may need more than
> MachinePrecision digits of working precision to meet these tolerances."
>
> I can see (by monitoring the Norm of the residual) that the algorithm
> gets trapped a long way from the minimum -- but I don't see why more
> digits of working precision are required.
>
> As an aside, is there an easy way (a handle?) to monitor the norm of the
> residuals at each step? One way is to write a function to compute the
> residuals,
>
> r[a_, b_, c_, d_] =
> Norm[(f[a, b, c, d] /@ data[[All,1]]) - data[[All,2]], Infinity];
>
> and then track them using StepMonitor or EvaluationMonitor
>
> Reap[FindFit[data, f[a, b, c, d][m], {a,b,c,d}, m,
> NormFunction :> (Norm[#,Infinity]&),
> EvaluationMonitor :> Sow[r[a,b,c,d]]]]
>
> but is there a _direct_ way of accessing the residuals which are
> computed by FindFit anyway?
>
> Cheers,
> Paul
>
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}}}
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].
Maxim Rytin
m.r at inbox.ru