Re: FindFit power law problem
- To: mathgroup at smc.vnet.net
- Subject: [mg117105] Re: FindFit power law problem
- From: Gary Wardall <gwardall at gmail.com>
- Date: Thu, 10 Mar 2011 06:06:49 -0500 (EST)
- References: <il50uv$nck$1@smc.vnet.net>
On Mar 8, 4:38 am, Bill Rowe <readn... at sbcglobal.net> wrote:
> On 3/7/11 at 5:48 AM, jgchild... 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.
Greg,
I was wondering if you really meant to transform both the x and y
data?
Fit[Log[10, data], {1, x}, x]
Does just that. Some software transforms just the y data when
estimating the coefficients for the model y = a x^b.
What happens is they transform just the y data and use the linear
model
yy=aa+bb*ln[x]
Finding the coefficients aa and bb and state the results as
y=Exp[aa]*x^bb
Since:
y=Exp[aa+bb*ln[x]]
is equivalent to:
y=Exp[aa]*Exp[bb*ln[x]]
in turn is equivalent to:
y=Exp[aa]*Exp[ln[x^bb]]
which is also equivalent to:
y=Exp[aa]*x^bb
Note then that a=Exp[aa] and b=bb.
I hope my Algebra is correct.
Gary