# Precision of LinearModelFit with Polynomials

I have a Problem regarding the fit of given points with a polynomial up to the fifth degree.

 tableofvalues=Import["tableofvalues.csv"]


My polynomial is:

polynomtabelle={1, x[1], x[1]^2, x[1]^3, x[1]^4, x[1]^5, x[2], x[2]^2, x[2]^3, x[2]^4, x[2]^5, x[3], x[3]^2, x[3]^3, x[3]^4, x[3]^5, x[4], x[4]^2,  x[4]^3, x[4]^4, x[4]^5}

polynom = LinearModelFit[tableofvalues, polynomtabelle,  Array[x, {4}]];


I know that my points have a certain symmetry. So the linear factors of coordinate x[1] have to be the same as for x[4]. The same applies for x[2] and x[3]. But if I ask for the coefficients I get for example:

NumberForm[Extract[CoefficientList[polynom["BestFit"], Array[x, {4}]], {5, 0, 0, 0} + 1], Floor@$MachinePrecision] 0.0619247206844402 NumberForm[Extract[CoefficientList[polynom["BestFit"], Array[x, {4}]], {0, 0, 0, 5} + 1], Floor@$MachinePrecision]
0.0619247173030426


One sees that they differ in the eigth digit after the decimal point, which is an accuracy too low for later treatment of the polynomial. Especially we have an old fitting routine written in FORTRAN in our working group that achieves 0.0619247173 for both coefficients. (Infact I intended to replace it with this one written in Mathematica).

How can I achieve numerically a better result. Perhaps you will ask for the precision of my input.

Precision[tableofvalues]
MachinePrecision


I run Mathematica on a 64Bit Computer so

$MachinePrecision 15.9546  The tableofvalues.csv is: 0.,0.,0.,0.,-0.9084906825 -0.04724315325,0.,0.,0.,-0.9075435543 -0.0377945226,0.,0.,0.,-0.9078643076 -0.028345891949999997,0.,0.,0.,-0.9081173929 -0.0188972613,0.,0.,0.,-0.9083049166 -0.00944863065,0.,0.,0.,-0.9084286871 0.00944863065,0.,0.,0.,-0.908492608 0.0188972613,0.,0.,0.,-0.9084363935 0.028345891949999997,0.,0.,0.,-0.9083237119 0.0377945226,0.,0.,0.,-0.9081562202 0.04724315325,0.,0.,0.,-0.9079356124 0.,-0.04724315325,0.,0.,-0.9074790503 0.,-0.0377945226,0.,0.,-0.9078199288 0.,-0.028345891949999997,0.,0.,-0.9080895811 0.,-0.0188972613,0.,0.,-0.9082900218 0.,-0.00944863065,0.,0.,-0.9084230317 0.,0.00944863065,0.,0.,-0.9084947071 0.,0.0188972613,0.,0.,-0.9084370745 0.,0.028345891949999997,0.,0.,-0.9083194184 0.,0.0377945226,0.,0.,-0.9081434864 0.,0.04724315325,0.,0.,-0.9079110427 0.,0.,-0.04724315325,0.,-0.9074790503 0.,0.,-0.0377945226,0.,-0.9078199288 0.,0.,-0.028345891949999997,0.,-0.9080895811 0.,0.,-0.0188972613,0.,-0.9082900218 0.,0.,-0.00944863065,0.,-0.9084230317 0.,0.,0.00944863065,0.,-0.9084947071 0.,0.,0.0188972613,0.,-0.9084370745 0.,0.,0.028345891949999997,0.,-0.9083194184 0.,0.,0.0377945226,0.,-0.9081434864 0.,0.,0.04724315325,0.,-0.9079110427 0.,0.,0.,-0.04724315325,-0.9075435543 0.,0.,0.,-0.0377945226,-0.9078643076 0.,0.,0.,-0.028345891949999997,-0.9081173929 0.,0.,0.,-0.0188972613,-0.9083049166 0.,0.,0.,-0.00944863065,-0.9084286871 0.,0.,0.,0.00944863065,-0.908492608 0.,0.,0.,0.0188972613,-0.9084363935 0.,0.,0.,0.028345891949999997,-0.9083237119 0.,0.,0.,0.0377945226,-0.9081562202 0.,0.,0.,0.04724315325,-0.9079356124 • Just curious: Adding degrees of freedom artificially will make your results worse. If you already know that the model's coefficients have a certain relationship, why don't you state that right in the model? – Dr. belisarius Jun 10 '15 at 18:36 • I made this table as minimal example. But I also need coupling constants lateron. For example the linear factor a*x[1]*x[4]. So it's not possible to just skipt x[3] and x[4] and mirror x[1] and x[2] onto them. For the diagonal elements this was of course possible. – mcocdawc Jun 10 '15 at 19:08 • Beside from this I would like to understand where this decrease of precision of 8 orders of magnitude comes from. – mcocdawc Jun 10 '15 at 19:09 • Using your data, the results from your NumberForm[ ] in my machine are more like 6.39247197093749 .... – Dr. belisarius Jun 10 '15 at 19:15 • Unfortunately I can't help you in pinpointing the origin of the drop in precision, but it seems that you can get what you want by manually setting your input values to an arbitrary (high) precision: tableofvalues = SetPrecision[#, 40]& @ {yourdata}. You will then obtain coefficients from your expression that are numerically the same to within the precision you have set. Further note: I also get the 6.39... answers from your code that @belisarius reported. – MarcoB Jun 10 '15 at 19:26 ## 1 Answer UPDATE A more accurate explanation than the culprit being "low variability" is that because all of the dependent variable values begin with "-0.907" or "-0.908" essentially "eats up"/"wastes" the first 3 significant digits. Simply subtracting the minimimum value of the dependent variable works even better than standardizing by the mean and standard deviation: (* Subtract minimum value from the dependent variable *) tableofvalues[[All, 5]] = tableofvalues[[All, 5]] - Min[tableofvalues[[All, 5]]];  The resulting values of the two coefficients in question become 0.0619247172727312 0.0619247172719015  This is a very common issue that can produce odd results for nearly all statistical packages. Looking intensely at the data and using some form of standardization is always recommended. End of UPDATE The main culprit is the low variability of the dependent variable: tableofvalues[[All,5]]. If you standardize that variable (subtracting the mean and then dividing by the standard deviation), then all your troubles go away. In fact standardizing variables is almost always recommended even for the independent variables whenever there might be round-off error involved especially for the fitting of polynomials. However, the suggestion by @MarcoB (setting a higher precision) achieves essentially the same result. stdvalues = tableofvalues; stdvalues[[All, 5]] = Standardize[tableofvalues[[All, 5]]]; stdev5 = StandardDeviation[tableofvalues[[All, 5]]]; polynomtabelle = {1, x[1], x[1]^2, x[1]^3, x[1]^4, x[1]^5, x[2], x[2]^2, x[2]^3, x[2]^4, x[2]^5, x[3], x[3]^2, x[3]^3, x[3]^4, x[3]^5, x[4], x[4]^2, x[4]^3, x[4]^4, x[4]^5} polynom = LinearModelFit[stdvalues, polynomtabelle, Array[x, {4}]]; NumberForm[ stdev5 Extract[ CoefficientList[polynom["BestFit"], Array[x, {4}]], {5, 0, 0, 0} + 1], Floor@$MachinePrecision]

NumberForm[
stdev5 Extract[
CoefficientList[polynom["BestFit"], Array[x, {4}]], {0, 0, 0, 5} +
1], Floor@$MachinePrecision]  with output (* 0.0619247172791471 *) (* 0.0619247172736994 *)  Note that one must multiply by the standard deviation of tableofvalues[[All,5]] to get the appropriate coefficient. • Thank you very much! If I apply the solution of @MarcoB, I got some interesting behaviour. With tableofvalues = SetPrecision[#,Floor@ MachinePrecision] &@Import["tableofvalues.csv"]; everything works fine BUT if I use SetPrecision[#,MachinePrecision] the problem with different values of symmetrical terms arises. – mcocdawc Jun 10 '15 at 23:17 • If I replace the top 3 lines of the standardization code above with stdvalues = Flatten[SetPrecision[#, 40] &@{tableofvalues}, 1];, then I get (after removing the multiplication by stdev5) both results being 0.0619247172725849. So that is likely a better result than doing the standardization that I suggested. – JimB Jun 11 '15 at 0:12 • I'd also like to suggest that fitting 21 coefficients with 41 data points, seems a bit excessive not to mention the use of a 5th degree polynomial for all independent variables. And from your comments it sounds like you're about to include additional interaction terms. You might consider using something like AIC to decide when overfitting becomes a problem. – JimB Jun 11 '15 at 0:22 • I wanted to avoid SetPrecision[#, 40] &@{tableofvalues} or any other arbitrary precision number to prevent numerically inefficient code. And the precision using Floor@$MachinePrecision is sufficient now. So my actual problem is solved. BUT again, the behaviour of SetPrecision is really strange. I get good results for every Integer number and hence also for Floor@MachinePrecision and Round@MachinePrecision. But as soon as I use SetPrecision[@,MachinePrecision] symmetrical terms differ in the eigth digit after decimal point. – mcocdawc Jun 11 '15 at 8:50
• This can be reproduced with the data given. @Jim Baldwin the tableofvalues posted is just a minimal example. Infact I use 361 points and 73 degrees of freedom. – mcocdawc Jun 11 '15 at 9:09