MathGroup Archive 2011

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

Search the Archive

Re: FindFit power law problem

  • To: mathgroup at
  • Subject: [mg117043] Re: FindFit power law problem
  • From: Bill Rowe <readnews at>
  • Date: Tue, 8 Mar 2011 05:37:14 -0500 (EST)

On 3/7/11 at 5:48 AM, jgchilders at (Greg Childers) wrote:

>I'm having a problem with FindFit and a power law problem.  Here's
>the data:

>data = {{1004, 0.003977}, {9970, 0.006494}, {100000, 0.012921},
>{1001000, .059795}}

>I'm wanting to fit it to a function of the form y = a x^b, and
>determine the best value of b.  When entered into Excel, it returns
>the exponent b = 0.383.  However, Mathematica gives

>FindFit[data, a x^b, {a, b}, x] {a->0.0000145749, b->0.601807}

>A graph of these values overlaid on the original data simply didn't
>look right.  Another way to find the exponent b is to take the log
>of both sides and do a linear fit:

>Fit[Log[10, data], {1, x}, x] -3.64936 + 0.383174 x

>And sure enough the exponent is 0.383 in agreement with Excel.  Why
>does FindFit give a different value?

The short answer is FindFit is trying to *correctly* solve the
following problem,

y = a x^b + error

where the error term is assumed to be normally distributed.
Excel most definitely is not doing this. Excel does precisely
what you did by taking the Log of your data then estimating the
parameters. Note, when you take the Log you have

Log[y] = Log[a] + Log[x^b + error]


Log[y] = Log[a] + b Log[x] + error

It is this last problem Excel solves and what FindFit solves
when you first take the Log of your data array then use FindFit.

An additional problem is seen when you do ListLogLogPlot[data].
That is, the resulting plot shows curvature indicating y = a x^b
isn't the underlying model for your data.

The curvature suggests a model like y = a(x+b)^c would better
fit your data. When I use this model I get

In[6]:= FindFit[data, a*(x + b)^c, {a, b, c}, x]

Out[6]= {a -> 2.14496457174612*^-6, b -> 32107.42765117816,
    c -> 0.7391222616795872}

and on checking

In[7]:= model = a (x + b)^c /. %

Out[7]= 2.14496*10^-6 (x+32107.4)^0.739122

The plot that results from

LogLogPlot[model, {x, 1000, 1.5 10^6},
  Epilog -> {PointSize[.01], Point[Log[data]]}]

seems to be in good agreement with the data

Note, I am not claiming this is a correct model for your data.
It is merely the first model that came to mind after looking a a
log log plot of your data points. And since you only have 4 data
points and this model has three free parameters, it is virtually
guaranteed to show a better fit to your data regardless of
whether it matches the "true" model or not. I've put true in
quotes here since if this data is experimental data, the "true"
model is often unknown.

  • Prev by Date: Re: FindFit power law problem
  • Next by Date: Re: evaluation-- one or many levels, your thoughts?
  • Previous by thread: Re: FindFit power law problem
  • Next by thread: Re: FindFit power law problem