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.