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