> restart:
> f := proc(x) 2*x + 1; end:
> rando := rand(-10^5 .. 10^5):
> for i from 1 to 10 do
>     x[i] := i;
>     y[i] := f(i) + rando()/(10.0^5)/2;
> od;
> 

                              x[1] := 1


                         y[1] := 3.159970000


                              x[2] := 2


                         y[2] := 4.938625000


                              x[3] := 3


                         y[3] := 7.277705000


                              x[4] := 4


                         y[4] := 9.361475000


                              x[5] := 5


                         y[5] := 11.13348500


                              x[6] := 6


                         y[6] := 12.98394000


                              x[7] := 7


                         y[7] := 14.80887500


                              x[8] := 8


                         y[8] := 17.03458000


                              x[9] := 9


                         y[9] := 19.46204500


                             x[10] := 10


                         y[10] := 21.04764000

> points := []:for i from 1 to 10 do
>     points := [op(points), [x[i], y[i]]];
> od:
> 
> points;

  [[1, 3.159970000], [2, 4.938625000], [3, 7.277705000],

        [4, 9.361475000], [5, 11.13348500], [6, 12.98394000],

        [7, 14.80887500], [8, 17.03458000], [9, 19.46204500],

        [10, 21.04764000]]

> p1 := plot(points, style=point,color = blue):
> A := Matrix(10,2): b := Vector(10):
> x_coords := []: y_coords:=[]: yw_coords:=[]:
> for i from 1 to 10 do
>     x_coords := [op(x_coords), x[i]]; A[i,1] := 1; A[i,2] := x[i];
>     y_coords := [op(y_coords), y[i]]; b[i] := y[i];
>     od:
> 
> A;

                              [1     1]
                              [       ]
                              [1     2]
                              [       ]
                              [1     3]
                              [       ]
                              [1     4]
                              [       ]
                              [1     5]
                              [       ]
                              [1     6]
                              [       ]
                              [1     7]
                              [       ]
                              [1     8]
                              [       ]
                              [1     9]
                              [       ]
                              [1    10]

> b;

                            [3.159970000]
                            [           ]
                            [4.938625000]
                            [           ]
                            [7.277705000]
                            [           ]
                            [9.361475000]
                            [           ]
                            [11.13348500]
                            [           ]
                            [12.98394000]
                            [           ]
                            [14.80887500]
                            [           ]
                            [17.03458000]
                            [           ]
                            [19.46204500]
                            [           ]
                            [21.04764000]

> x_coords;

                   [1, 2, 3, 4, 5, 6, 7, 8, 9, 10]

> y_coords;

  [3.159970000, 4.938625000, 7.277705000, 9.361475000, 11.13348500,

        12.98394000, 14.80887500, 17.03458000, 19.46204500,

        21.04764000]

> model := stats[fit,leastsquare[[x,y],y='c'[1]*x + 'c'[0]]]([x_coords,
> y_coords]);

               model := y = 1.997757576 x + 1.133167333

> with(LinearAlgebra):
> xhat := LinearSolve(Transpose(A) . A, Transpose(A) . b);

                            [1.13316733333333963]
                    xhat := [                   ]
                            [1.99775757575757496]

> Norm(A . xhat - b, 2);

                         0.601892010750393114

> p2 := plot(op(2, model), x=0..10, color = red):
> plots[display]({p1, p2});

>