Re: FindFit power law problem
- To: mathgroup at smc.vnet.net
- Subject: [mg117043] Re: FindFit power law problem
- From: Bill Rowe <readnews at sbcglobal.net>
- Date: Tue, 8 Mar 2011 05:37:14 -0500 (EST)
On 3/7/11 at 5:48 AM, jgchilders at mailaps.org (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]
not
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.