How to use LeastSquares[] or QRDecomposition[] to Fit y = mx + b

Questions
, ,
1

I am trying to calculate the m and b values as a best fit to the equation of the line of y=mx + b - see attached notebook.

In the first section of the notebook, I performed the detail calculations manually to obtain the variables of m = 1.1 and b = 1 to create the equation of the line of y = 1.1x + 1.

Then I tried to use two Mathematica functions (QRDecomposition[] and LeastSquares[]) in an attempt to streamline the calculations and ran into problems. I tried to use examples in both Mathematica as well as examples in the WolframU course on Linear Algebra in which I managed to thoroughly confuse myself because none of the new calculations even remotely resembles the calculations or the end results that I performed manually.

If there is one, I would certainly appreciate someone showing me the recommended method of using Mathematica functions to perform the calculations that I performed manually in the first section of the attached notebook.

Thank you,

Mitch Sandlin

QR_Fact.nb (23.2 KB)

2

Here’s how I would do it. (I made up my own data vecs xx, yy.)

xx = Subdivide[0., 5., 10];
yy = 1 + 1.1*xx + 
   RandomVariate[NormalDistribution[0., 0.2], Length@xx];
design = Transpose@{ConstantArray[1., Length@xx], xx};
bmvec = LeastSquares[design, yy]
(*  {0.936894, 1.11642}  *)

{qq, rr} = QRDecomposition[design];
LinearSolve[rr, qq . yy]
(*  {0.936894, 1.11642}  *)

ListPlot[Transpose[{xx, yy}], 
 Epilog -> {Red, Line[{{0, First[bmvec]}, {5, bmvec . {1, 5}}}]}] (*Edit: was Last[bmvec]*)

1882Untitled.png
1882Untitled.png1080×696 19.5 KB

Edit note: In editing the code for posting, I switched the order from {m, b} to {b, m} in what is now called bmvec; but I neglected to change Last to First in ListPlot[].

3

You can do it with Minimize or FindFit:

y = {2, 3, 5, 5};
x = {{1, 1}, {2, 1}, {3, 1}, {4, 1}};
\[Beta] = {m, b};
Minimize[(y - x . \[Beta]) . (y - x . \[Beta]), {m, b}]
data = Transpose[{x[[All, 1]], y}];
Block[{x}, FindFit[data, m*x + b, {m, b}, x]]
4

Hi Michael;

I find your solution very interesting. In studying your solution, I tried to make modifications to incorporate the two matrices that I used in my initial notebook, with some success. The calculations seem to work fine, but the ListPlot[] function fails. I have tried several modifications without any success and am not sure what is causing my problem. - please see attached.

Thanks,

Mitch Sandlin

QR_Fact3x.nb (11.6 KB)

5

This gives a result:

ListPlot[Transpose[{xx[[All, 1]],
   yy[[All, 1]]}],
 Epilog -> {Red,
   Line[{Flatten@{0, Last[bmvec]},
     {5, Flatten@bmvec . {1, 5}}}]}]
6

Some differences in setup:

xx
(* 
{{1, 1}, {2, 1}, {3, 1}, {4, 1}}  <-- Yours is a design matrix (rank 2 array)
{1, 2, 3, 4}                      <-- Mine is a vector of x abscissae (rank 1 array)
*)

yy
(*
{{2}, {3}, {5}, {5}}  <-- Yours is a column matrix
{2, 3, 5, 5}          <-- Mine is a vector
*)

If you want to set up the problem all in matrices, then the solution will be in terms of matrices. The plot can be adapted to the matrices as follows:

ListPlot[Transpose[{xx[[All, 1]], yy[[All, 1]]}], 
 Epilog -> {Red, 
   Line[{{0, First@First@bmvec}, (* there was a typo in my previous code *)
    {5, {{1, 5}} . bmvec // First // First}}]}]

The code can be fixed in many ways. The following couplets are pairwise equivalent:

Transpose[{xx[[All, 1]], yy[[All, 1]]}]
Table[{xx[[1 k, 1]], yy[[k, 1]]}, {k, Length@xx}]

{{1, 5}} . bmvec // First // First
({{1, 5}} . bmvec)[[1, 1]]

First@First@bmvec
bmvec[[1, 1]]

I’ve found it easier to work with vectors as rank-1 arrays and not as rank-2 ones (matrices) in Mathematica. For one thing, it’s how the examples in the documentation are set up.

†The rank of an array corresponds to how many layers of curly braces the expression has. It also corresponds to how many indices are needed between the square brackets [[…]] to extract an entry.

7

Hi Gianluca;

In the Minimize function that you suggested, you used vector representation to describe the data. I intern tried to perform the same calculations except using Matrice representation for the data and, in spite of my best efforts, could not get the Minimize function to calculate. I tried using the functions: Flatten[], Transpose[], and other methods to get the calculations to work all without success for either the Minimize[] or ListPlot[] functions - please see attached notebook.

As an aside, I cannot understand why simply representing data in a different format makes the calculations and data display so different and I hope that someone can explain what programming changes need to be made. On the surface, it appears that the vector form may be the less complicated method to process data, However, it is my understanding that in Tensor Calculus, the data is represented as Matrices. In any event, I am thoroughly confused and hope you can shed some light,

Thanks so much,

Mitch Sandlin

QR_FactMin.nb (43.7 KB)

8

If the objective is to find a fit, I suggest using LinearModelFit. Why? (1) It’s simpler than what you’ve tried. (2) You get measures of precision (standard errors) automatically. (3) You get access to various fit diagnostic statistics.

data = {{1, 2}, {2, 3}, {3, 5}, {4, 5}}
lmf = LinearModelFit[data, x, x];
lmf["ParameterTable"]

parametertable.png

Show[ListPlot[data], Plot[lmf[x], {x, 1, 4}]]

figure.png
figure.png528×340 11.3 KB

9

Hi Glanluca, Michael & Jim;

I want to thank all of you for the assistance you gave me. With your help, I finally figured out the solution, and the journey was certainly no walk in the park. However, now I have a much better understanding of the differences between matrices and vectors and now understand that a little different skill set is needed when working with matrices.

Thanks again and have a great day,

Mitch Sandlin