# LinearModelFit with Standardized data

I have a small data set of six points:

data1={{2014.,0.015},{2015.,0.005},{2016.,0.0},{2017.,0.01},{2018.,0.02},{2019.,0.014}};

ListPlot[data1
,Frame->True
,PlotRange->{{2013.,2022.},{-0.01,0.03}}
,PlotStyle->Directive[Orange,PointSize[Large]]
]


giving the following plot

The Mean and Standard Deviation of data1 are

Mean[data1]
StandardDeviation[data1]


{2016.5,0.0106667}

{1.87083, 0.00725718}

The data1 can be transformed to data2, having zero Mean and unit Standard Deviation, by utilizing Standardize[.]

data2=Standardize[data1];

ListPlot[data2
,Frame->True
,PlotRange->{{-2.,2.},{-2.,2.}}
,PlotStyle->Directive[Orange,PointSize[Large]]
]


In any case, LinearModelFit[.] allows to fit a polynomial trough the data

lmFit[data_List,degree_Integer]:=LinearModelFit[data,Table[x^i,{i,degree}],x]


Mathematically, a 5th degree polynomial fits exactly through any six data points. However, trying to fit a polynomial of degree=5 to gives quite different results

lmFitPlot[data_List,degree_Integer,{xmin_,xmax_,ymin_,ymax_}]:=Module[{lmf,ss},

lmf=lmFit[data,degree];
ss=Total[lmf["FitResiduals"]^2]; (* Sum of squared residuals *)

Show[
{Plot[lmf[x],{x,xmin,xmax}
,PlotRange->{{xmin,xmax},{ymin,ymax}}]
,ListPlot[data,PlotStyle->Directive[Orange,PointSize[Large]]]
}
,Frame->True
,FrameLabel->{{"",""},{"Year",Row[{"Sum squared residuals= ",ss}]}}
,ImageSize->Medium
]
]


For data1

lmFitPlot[data1,5,{2013.,2022.,-0.01,0.03}]


gives a very noisy fit

While for data2 the fit is quite decent

lmFitPlot[data2,5,{-2.,2.,-2.,2.}]


This discrepancy is caused by the poor rank of the Design Matrix for data1

MatrixRank[lmFit[data1,5]["DesignMatrix"]]
MatrixRank[lmFit[data2,5]["DesignMatrix"]]


3

6

Since there exists is a straightforward Geometric Transformation between data1 and data2

FindGeometricTransform[data1,data2]


my question is: Does there exist an inverse geometric transformation which transforms the decent linear regression model of data2 back to the original coordinate system of data1?

Wikipedia shows a few interesting transformation examples Geometric transformation

Thanks.

• Knowing the TransformationFunction you could use InverseFunction (see documantation ) – Ulrich Neumann Apr 23 '19 at 10:22
• Ulrich Neumann, Thanks for your reply. Yes, I can make the InverseFunction of the Geometric transformation. But how do I transform the polynomial function from the LinearModelFit? I could only find transformations of geometric objects: circle, polygon, etc. – Romke Bontekoe Apr 23 '19 at 14:10
• @ Romke Bontekoe The transformation matrix is a blockmatrix of the form {{A,b},{c^t,1}}. The rational transformation point x->point y is defined as y=(A.x+b)/(c.x+1) – Ulrich Neumann Apr 23 '19 at 14:15
• @Ulrich Neumann, I would much appreciate if you would take the lmFit[data2,5] solution (from the above code) and transform this polynomial as you describe. I failed to do so myself. Please post this as an answer. – Romke Bontekoe Apr 24 '19 at 6:14

As appreciated here my answer concerning the GeometricTransformation

First let's define

 y2 = Normal[lmFit[data2, 5]] /. x -> x2(* data2 coordinate system {x2,y2}*)


and the geometric transformation

 gt12 = FindGeometricTransform[data1, data2][[2]](* data1~gt12[data2]*)


which transforms {x2,y2} to {x1,y1}(coordinate system data1).

Now the final transformation back to the {x1,y1} coordinates:

{x1, y1} = gt12[{x2, y2}] // Simplify
Show[{ParametricPlot[{x1, y1}, {x2, -1.33631, 1.33631}],Graphics[{Red, Point[data1]}]}, AspectRatio -> 1]


which fit's the data.

Knowing gt12 you can also describe the inverse transformation

gt21=InverseFuncion[gt12] (* gt21[data1]~data2 *)
Chop[ gt21[data1] - data2, 10^-9 ]
(*{{0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}, {0, 0}}*)

• Ulrich Neumann, thank you for this insightful solution. I have learned a lot. – Romke Bontekoe Apr 24 '19 at 15:58
• You're welcome! – Ulrich Neumann Apr 24 '19 at 19:44
• I have one small question still. How did you derive the ParametricPlot interval {x2, -1.33631, 1.33631} from the data endpoints? I am still trying to understand the mechanism of the geometric transformations. – Romke Bontekoe Apr 25 '19 at 6:01
• These values I took from data2 data range! – Ulrich Neumann Apr 25 '19 at 6:12
• I can compute gt21[{2014,0}] and gt21[{2019,0}] and take #[[1]] of these. But maybe you have a more elegant way. – Romke Bontekoe Apr 25 '19 at 6:44

You really have no business fitting or predicting a 5th degree polynomial with 6 data points. To add insult to injury your data has 2 sets of 3 points that are almost perfectly colinear. Finally, you'll have no estimate of error. Hopefully, this is just for a class exercise.

The anomalies that you see are due to a lack of precision. That can be remedied by rationalizing your data and increasing the WorkingPrecision for LinearModelFit and Plot.

data1 = {{2014., 0.015}, {2015., 0.005}, {2016., 0.0}, {2017., 0.01}, {2018., 0.02}, {2019., 0.014}};
data1 = Rationalize[data1, 0];

lmFit[data_List, degree_Integer] := LinearModelFit[data, Table[x^i, {i, degree}], x,
WorkingPrecision -> 50]

lmFitPlot[data_List, degree_Integer, {xmin_, xmax_, ymin_, ymax_}] :=
Module[{lmf, ss}, lmf = lmFit[data, degree];
ss = Total[lmf["FitResiduals"]^2];(*Sum of squared residuals*)
Show[{Plot[lmf[x], {x, xmin, xmax},
PlotRange -> {{xmin, xmax}, {ymin, ymax}},
WorkingPrecision -> 50],
ListPlot[data, PlotStyle -> Directive[Orange, PointSize[Large]]]},
Frame -> True,
FrameLabel -> {{"", ""}, {"Year",
Row[{"Sum squared residuals= ", ss}]}}, ImageSize -> Medium]]

lmFitPlot[data1, 5, {2013, 2022, -0.01, 0.03}]


There's no need to standardize and then transform back in this case:

data2 = Standardize[data1];
lmFitPlot[data2, 5, {-2., 2.5, -2., 2.}]


In this case if you don't take care of the precision issues in the first place, transforming back from the standardized model you might still have the precision issues.

Probably the simplest approach is just to subtract 2014 from the "x" values. Then there's no need to add in a WorkingPrecision statement.

data1 = {{2014., 0.015}, {2015., 0.005}, {2016., 0.0}, {2017., 0.01}, {2018., 0.02}, {2019., 0.014}};
data3 = data1;
data3[[All, 1]] = data3[[All, 1]] - 2014;

lm3 = LinearModelFit[data3, Table[x^i, {i, 5}], x];
Show[Plot[lm3[x - 2014], {x, 2013.5, 2021}],
ListPlot[data1, PlotStyle -> {Red, PointSize[0.02]}]]


• +1 simply for that opening sentence. I found that rather amusing. – Q.P. Apr 23 '19 at 15:31
• @JimB, thanks for your answer. You guessed right, I am developing this example for educational purposes. Most likely, the targeted audience does not have access to arbitrarily high working precisions. Therefore your first solution is not in the right direction, though it is very instructive to me. Your additional solution is exactly overlapping your first solution, to within 10^-15, which I did not expect from an arbitrary shift of the x-values by subtracting 2014. It might be interesting to you that these are real financial ratio data. Such data never come with error bars. Thanks. +1 – Romke Bontekoe Apr 24 '19 at 8:03
• "Such data never come with error bars." Maybe I'm misunderstanding but it's rare that data comes with error bars (from any scientific field). However, one can still estimate the "error" as the variance about the curve is one of the parameters estimated in the fitting of a regression (except when one fits a 5-th degree polynomial with 6 data points. And NOT stating estimates of error when one can is unfortunately not limited to economists/financial folks. – JimB Apr 24 '19 at 14:55
• The lack of precision when using just the raw data comes from, for example, raising 2019 to the 5th power: 2019^5 = 33,549,155,665,686,099. Necessary digits get lost very quickly. After subtracting 2014 from the predictors the largest value is (2019-2014)^5 = 3125 and there is little (if any) loss of precision. Because economists deal with current years and high order polynomials, maybe this is the source of the joke: "Economists were invented to make the meteorologists feel better about their predictions." – JimB Apr 24 '19 at 15:05
• JimB, my plan is to demonstrate the numerical instabilities, which many professionals are unaware of (unfortunately). I want to treat the unknown uncertainty in the data by employing a so called Jeffreys' prior, as part of a Bayesian treatment of this regression problem. Thanks once more. – Romke Bontekoe Apr 24 '19 at 16:05