Re: FindFit power law problem
- To: mathgroup at smc.vnet.net
- Subject: [mg117156] Re: FindFit power law problem
- From: Ray Koopman <koopman at sfu.ca>
- Date: Thu, 10 Mar 2011 06:16:03 -0500 (EST)
- References: <il2d8k$65a$1@smc.vnet.net> <il8r0t$lvt$1@smc.vnet.net>
On Mar 9, 1:21 pm, Peter Pein <pet... at dordos.net> wrote:
> Am 07.03.2011 11:49, schrieb Greg Childers:
>> Hi,
>>
>> 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?
>>
>> Greg
>
> Hi Greg,
>
> if you really insist on b being 0.383, then find the weights:
>
> In[1]:=
> data={{1004,0.003977},{9970,0.006494},{100000,0.012921},{1001000,.059795}};
>
> In[2]:=
> bw[w_List/;(And@@NumericQ/@w)]:=Cases[Normal[NonlinearModelFit[data,a
> x^b,{a,b},x,Weights->w]],Power[x,b_]:>b,1,1][[1]]
> In[3]:= weights=w/@Range[4];
> In[4]:=
> soln=NMinimize[{(bw[weights]-0.383)^2,And@@Flatten[{Thread[weights>=0],Norm[weights]==1,weights[[1]]<=weights[[2]],weights[[3]]>=weights[[4]]
> (* put focus on the 'middle'*)}]},weights]
> Out[4]=
> {2.83106*10^-14,{w[1]->0.264427,w[2]->0.964382,w[3]->0.00681409,w[4]->0.000760789}}
> In[5]:= NonlinearModelFit[data,a
> x^b,{a,b},x,Weights->((w/@Range[4])/.soln[[2]])]//Normal
> Out[5]= 0.000196237 x^0.383
>
> Peter
Here is the same analysis, but done by version 5.2. It gives the
desired 'b', but 'a' is only loosely the same, and the weight
distribution is flatter. This is a nice example of the extent to which
the weights in weighted least squares fitting are not always critical.
In[1]:= <<Statistics`NonlinearFit`
In[2]:= data = {{1004,0.003977},{9970,0.006494},{100000,0.012921},
{1001000,.059795}};
In[3]:= bw[w_List/;(And@@NumericQ/@w)] := Cases[Normal[
NonlinearFit[data, a x^b, x, {a,b}, Weights->w] ],
Power[x,b_] :> b, 1, 1 ][[1]]
In[4]:= weights = w /@ Range[4];
In[5]:= soln = NMinimize[{(bw[weights]-0.383)^2, And@@Flatten[
{Thread[weights>=0], Norm[weights]==1,
weights[[1]]<=weights[[2]], weights[[3]]>=weights[[4]]
(* put focus on the 'middle'*)}]}, weights]
Out[5]= {3.9472*^-22, {w[1] -> 0.590995, w[2] -> 0.806552,
w[3] -> 0.0140113, w[4] -> 0.00176326}}
In[6]:= NonlinearFit[data, a x^b, x, {a,b},
Weights->((w/@Range[4])/.soln[[2]])]//Normal
Out[6]= 0.000203453 x^0.383