Re: Strange results for simple calculations
- To: mathgroup at smc.vnet.net
- Subject: [mg109535] Re: Strange results for simple calculations
- From: Bill Rowe <readnews at sbcglobal.net>
- Date: Tue, 4 May 2010 06:29:08 -0400 (EDT)
On 5/3/10 at 6:10 AM, astronerma at gmail.com (koringkriek) wrote:
>Hi All, I have a data which I fit with the polynomial of 4th order
>and then I want to find the mean of the differences between the data
>and the model:
>(* my data, 21 points *) data = {{47.5, 0.01484}, {48.5, 0.01603},
>{49.5, 0.017195}, {50.5, 0.0240}, {51.5, 0.033409}, {52.5, 0.0467},
>{53.5, 0.05748}, {54.5,0.063885}, {55.5, 0.071935}, {56.5, 0.0748},
>{57.5, 0.07639}, {58.5, 0.07639}, {59.5, 0.0748}, {60.5, 0.069635},
>{61.5, 0.0609}, {62.5, 0.05578}, {63.5, 0.05167}, {64.5, 0.0450},
>{65.5, 0.0406}, {66.5, 0.0381}, {67.5, 0.03641}};
>(* polynomial fit *) f = Fit[data, {1, x, x^2, x^3, x^4}, x]
This is not the best way to fit a 4th order polynomial to your
data. The problem with this approach is the powers of x do not
form an orthogonal set of basis functions. Additionally, using
either Fit or FindFit with powers of x in this fashion is an
increasingly unstable numeric problem as the order of the
polynomial increases. A far better choice would be to use
orthogonal polynomials such as Chebyshev polynomials.
>(* differences between the data an model *) var = Table[data[[i, 2]]
>- f /. x -> data[[i, 1]], {i, 1, Length[data]}]
I would compute this as:
var = data[[All, 2]] - (f /. x -> data[[All, 1]]);
or
var = (Last@# - (f /. x -> First@#)) & /@ data;
rather than using Table to access each element of data individually.
>(* calculate mean *) Mean[var]
>My result is: -3.65422*10^-14 which I find "a bit" underestimated.
With:
In[31]:= $Version
Out[31]= 7.0 for Mac OS X x86 (64-bit) (February 19, 2009)
I get:
In[32]:= Mean[var]
Out[32]= 1.24165*10^-13
The difference between my result and yours is certain to be the
result of using either different versions of Mathematica or
different platforms with machine precision arithmetic.
>Same goes wrong if I just want to do Sum[var]. Now compare the two
>loops where the difference is only in the upper limit of the sum 20
>or 21:
>(* loop 1 *) i = 1; a = 0; While[i <= 21,
>a = a + var[[i]]; i = i + 1;
>]
>a
>(* loop 2 *) i = 1; a = 0; While[i <= 20,
>a = a + var[[i]]; i = i + 1;
>]
>a
>First loop returns -7.67386*10^-13 the second 0.00251855. So, is
>there something wrong with the last point in my table var? How about
>I change the initial value of i and let it count till 21,
Loops really aren't needed. Your first loop sums all of the
elements of var while your second loop sums all but the last. Since,
In[33]:= Last@var
Out[33]= -0.00251855
and
In[34]:= Total[Most@var]
Out[34]= 0.00251855
it should not be surprising to get a very small number for the
sum of all elements of var.
>(* loop 3 *) i = 2; a = 0; While[i <= 21,
>a = a + var[[i]]; i = i + 1;
>]
>a
>And the result is: -0.000761973! I am completely lost here. Could
>someone please explain me what I am doing wrong?
This last loop sums all but the first element of var. So,
In[35]:= First@var
Out[35]= 0.000761973
In[36]:= Total[Rest@var]
Out[36]= -0.000761973
Again, the sum of all of the elements should be a small number.
In fact, since the sum of the negative elements is nearly equal
to the sum of the positive elements, deleting any one element
from the sum will give similar results.
>Could someone please explain me what I am doing wrong?
You really haven't done anything wrong. The explanation lies in
the values of the differences and the vagaries of machine
precision arithmetic. It isn't the loops you've used that cause
the apparent difficulty. However, the loops aren't necessary
aren't the most efficient way to use Mathematica.
Also, the reason for not using Fit in the manner you have done
above is the numerical instability that arises in this type of
problem for higher order polynomials. For myself, I would not do
polynomial fitting in this manner for any polynomial with powers
greater than 2.