> 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:
> 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.13316733333333586]
                    xhat := [                   ]
                            [1.99775757575757540]

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

                         .601892010750393114

# One can compute the infinity-norm nearest vector to b in the
# column-space of A:
> read("/afs/eos.ncsu.edu/users/k/kaltofen/www/courses/LinAlgebra/Maple/
> LeastSqus/incsolve.mpl"):
> xhat_inf := incsolve(convert(A,'matrix'), convert(b,vector));

                           w = .3064710000


                xhat_inf := [.9745480000, 2.020114000]


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

                         .349059484848485368

# A . xhat_inf is closer to b when using infinity norm:
> Norm(A . Transpose(convert(xhat_inf,'Vector')) - b, infinity);

                         .306471000000000160

> for i from 1 to 10 do
>     y_inacc[i] := f(i) + i^2*rando()/(10.0^5)/10;
> od;

                      y_inacc[1] := 2.953602000


                      y_inacc[2] := 5.147176000


                      y_inacc[3] := 7.041220000


                      y_inacc[4] := 10.57302400


                      y_inacc[5] := 11.73810000


                      y_inacc[6] := 10.53814000


                      y_inacc[7] := 19.57341500


                      y_inacc[8] := 20.24723200


                      y_inacc[9] := 25.08852700


                      y_inacc[10] := 25.42480000

> y_inacc_coords:=[]: y_weighted_coords:=[]:
> for i from 1 to 10 do
>     y_inacc_coords := [op(y_inacc_coords), y_inacc[i]];
>     y_weighted_coords := [op(y_weighted_coords), Weight(y[i],1/i)];
> od:
> model_inacc := stats[fit,leastsquare[[x,y],y='c'[1]*x +
> 'c'[2]]]([x_coords, y_inacc_coords]);

            model_inacc := y = 2.628251588 x - .6228601333

# The model is inaccurate because the errors grow at larger value of i. 
# A correction is possible by weighted norms:

> y_weighted_coords;

  [Weight(3.159970000, 1), Weight(4.938625000, 1/2),

        Weight(7.277705000, 1/3), Weight(9.361475000, 1/4),

        Weight(11.13348500, 1/5), Weight(12.98394000, 1/6),

        Weight(14.80887500, 1/7), Weight(17.03458000, 1/8),

        Weight(19.46204500, 1/9), Weight(21.04764000, 1/10)]

> model_weighted := stats[fit,leastsquare[[x,y],y='c'[1]*x +
> 'c'[2]]]([x_coords, y_weighted_coords]);
>               
> 

          model_weighted := y = 1.997472352 x + 1.134736066

> 
>