MathGroup Archive 2002

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

Search the Archive

Re: Problems with NonlinearRegress

  • To: mathgroup at smc.vnet.net
  • Subject: [mg32248] Re: Problems with NonlinearRegress
  • From: Adam Smith<adam.smith at hillsdale.edu>
  • Date: Wed, 9 Jan 2002 03:18:02 -0500 (EST)
  • References: <a1633l$pp0$1@smc.vnet.net>
  • Sender: owner-wri-mathgroup at wolfram.com

This was a very interesting study for me and the use of non-linear fits in
Mathematica. I apologize for the length of this, but I wanted to point out all
the pitfalls one can encounter in doing fits of this type.

Since I am not sure what your data "test4.txt" looked like I created a data set
based on your function as follows.


func = base + amplitude/(((freq - wo)/gamma)^2 + 1)

Note that I avoided the capital letters on "base" and "gamma".  Capital letters
are often reserved for functions and special values.

Then to create the data with an error of 1% in each point.

mysub = {base -> 0.65, gamma -> 0.008, amplitude -> 1.31, wo -> 8.82}
temprange = 0.008
data = Table[{freq, func*(1 + Random[Real, {-.01, 0.01}]) /. mysub}, {freq, 
8.82 - 15*temprange, 8.82 + 15*temprange, temprange/2.5}];
Clear[temprange];

Now check that the fake data looks OK

dataplot = 
ListPlot[data, PlotStyle -> PointSize[0.02], PlotRange -> All, 
DisplayFunction -> Identity];
exactplot = 
Plot[func /. mysub, {freq, 8.65, 8.95}, PlotRange -> All, 
DisplayFunction -> Identity];
Show[exactplot, dataplot, DisplayFunction -> $DisplayFunction];

The first problem that I see is that your function as defined is very sensitive
to the starting value of w0.  I was able to get reasonable results with your
other starting parameters only for starting values of w0 = 8.81 to 8.83


fit2 = NonlinearRegress[data, 
func, {freq}, {{base, 0.8}, {amplitude, 1}, {wo, 8.81}, {gamma, 0.009}}, 
MaxIterations -> 100, RegressionReport -> ParameterCITable] 

Estimate	Asymptotic SE
base		0.649568        0.000711214
amplitude	1.30892	        0.00353695
wo		8.81999	        0.0000216297
gamma		-0.00801306	0.0000342016

Anything less that 8.81 and greater than 8.83 produces unrealistic values of all
parameters with very large Asymptotic SE.  This makes sense if you look at the
distribution.  You need a starting value of wo that is at lest somewhere in the
"central peak".  Also note that the gamma parameter is negative but of the
correct magnitude.  This is because of the square on gamma in the function.

Now I found that a modification of the function definition produces a more
robust routine.  Based upon the definition I found of the Lorentzian:

1/Pi (g/2)/( (x-m)^2 + (g/2)^2)

I took you definition and massaged it a bit with your gamma = g/2 to get:

myfunc = base + amplitude*(gamma/((freq - wo)^2 + gamma^2))

Here I absorbed the 1/Pi into the "amplitude" so instead of 1.31 my exact
amplitude is 1.31*.008 = .01048.

Anyhow, I checked that this gave the same data as your function and it did.
Here the fit was much less sensitive to the starting parameters and worked for
starting values of w0 from 8.78 to 8.86 (vs. 8.81 to 8.83).  You still need to
be careful in the starting value, but now it seems to me this is a much easier
range to "eyeball" to pick the value.

myfit = NonlinearRegress[mydata, 
myfunc, {freq}, {{base, 0.8}, {amplitude, 1.*.009}, {wo, 8.8}, {gamma, 
0.009}}, MaxIterations -> 100] 
Estimate	Asymptotic SE
base		0.65063         0.00069442
amplitude	0.0105122       0.0000343314
wo		8.82            0.0000210092
gamma		0.00799454      0.000033212

Hope this helps.

In article <a1633l$pp0$1 at smc.vnet.net>, Martin Duemling says...
>
>Hi,
>I have the following problem: 
>I try to do a Lorentzian-fit with my data. It works very good if the datas are 
>simulated and very close to the ideal curve. But as soon as I use my real data 
>(they have still Lorentzian shape) I don't get any result which makes sense 
>(See example).
>
>Maybe somebody can help, Thank you.
>Martin
>
>
>Input:
>
>Clear["Global`*"]
>Needs["Statistics`NonlinearFit`"]
>data = ReadList["test4.txt", {Number, Number}]; 
>fit2 = NonlinearRegress[data, \
>    Base + Amplitude/((1 + ((((freq - wo))/gamma))^2)),
>    freq, {{Base, 0.8}, {Amplitude, 1}, {wo, 9}, {gamma, 0.009}},
>    Weights -> Equal, MaxIterations -> 100, WorkingPrecision -> 16,
>    RegressionReport
>       -> {ParameterTable,  ANOVATable, EstimatedVariance},
>    AccuracyGoal -> 16, PrecisionGoal -> 16, ShowProgress -> True]
>
>
>Output
>Iteration:1 ChiSquared:37.193  Parameters:{0.8, 1., 9., 0.009}
>
>Iteration:2 ChiSquared:37.5594  Parameters:{0.863029, 17.9184, 11.1608, 
>0.0648881}
>
>Iteration:3 ChiSquared:38.099  Parameters:{0.878984, 46.7423, 16.506, 
>0.117541}
>
>Iteration:4 ChiSquared:38.7406  Parameters:{0.900703, 29.929, 22.6742, 
>0.096431}
>
>Iteration:5 ChiSquared:38.7053  Parameters:{0.901547, -25.3515, 44.1125, 
>0.0073759}
>...
>And than the values are staying approx constant
>
>
>
>For comparison the correct values are (other fit program):
>Base= 0.65 
>Amplitude= 1.31
>wo= 8.82 
>gamma= 0.008
>
>
>

Adam Smith
Dept. of Physics
Hillsdale College
adam.smith at hillsdale.edu


  • Prev by Date: Re: How to call BinomialDistribution from within a package?[2]
  • Next by Date: Re: ORDINARY DIFFERENTIAL EQUATION
  • Previous by thread: Problems with NonlinearRegress
  • Next by thread: Wrestling with Mathematica on Trig Simplification