Precision of LinearModelFit with Polynomials

Question

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

Answer 1

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.