> f := proc(x) 2*x + 1; end;

                       f := proc(x) 2*x + 1 end

> rando := rand(-10^5 .. 10^5);

  rando := proc()
local t;
global _seed;
    _seed := irem(427419669081*_seed, 999999999989);
    t := _seed;
    irem(t, 200001) - 100000
end

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

                              x[1] := 1


                         y[1] := 2.953602000


                              x[2] := 2


                         y[2] := 5.036794000


                              x[3] := 3


                         y[3] := 7.004580000


                              x[4] := 4


                         y[4] := 9.098314000


                              x[5] := 5


                         y[5] := 11.02952400


                              x[6] := 6


                         y[6] := 12.93161500


                              x[7] := 7


                         y[7] := 15.09333500


                              x[8] := 8


                         y[8] := 17.05073800


                              x[9] := 9


                         y[9] := 19.07516700


                             x[10] := 10


                         y[10] := 21.04424800

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

  [[1, 2.953602000], [2, 5.036794000], [3, 7.004580000],

        [4, 9.098314000], [5, 11.02952400], [6, 12.93161500],

        [7, 15.09333500], [8, 17.05073800], [9, 19.07516700],

        [10, 21.04424800]]

> p1 := plot(points, style=point,color = blue):
> A := linalg[matrix](10,2); b := linalg[vector](10);

                   A := array(1 .. 10, 1 .. 2, [])


                       b := array(1 .. 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];
>     yw_coords := [op(yw_coords), Weight(y[i],1/i)]; b[i] := y[i];
> od:
> 

                            x_coords := []


                            y_coords := []


                           yw_coords := []

> print(A);

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

> print(b);

  [2.953602000, 5.036794000, 7.004580000, 9.098314000, 11.02952400,

        12.93161500, 15.09333500, 17.05073800, 19.07516700,

        21.04424800]

> x_coords;

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

> y_coords;

  [2.953602000, 5.036794000, 7.004580000, 9.098314000, 11.02952400,

        12.93161500, 15.09333500, 17.05073800, 19.07516700,

        21.04424800]

> model := stats[fit,leastsquare[[x,y],y=a*x + b]]([x_coords,
> y_coords]);

               model := y = 2.007287085 x + .9917127333

> with(linalg):
Warning, new definition for norm
Warning, new definition for trace
> linsolve(evalm(transpose(A) &* A), evalm(transpose(A) &* b));

                      [.9917127252, 2.007287086]

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

> read("/afs/eos.ncsu.edu/users/k/kaltofen/www/courses/LinAlgebra/Maple/
> LeastSqus/incsolve.mpl");

  incsolve := proc(A, b)
local m, n, i, j, c, d, w, x, s, ss;
    if type(A, matrix) = false then
        ERROR(`first parameter must be of type matrix`)
    fi;
    if type(b, vector) = false then
        ERROR(`second parameter must be of type vector`)
    fi;
    m := linalg[rowdim](A);
    n := linalg[coldim](A);
    if m < n then ERROR(`the number of rows must be > or =\
         the number of columns`)
    fi;
    if m <> linalg[vectdim](b) then
        ERROR(`incompatible dimensions`)
    fi;
    x := linalg[vector](n);
    d := linalg[vector](m);
    c := {};
    for i to m do
        d[i] := -b[i];
        for j to n do d[i] := d[i] + A[i, j]*x[j] od;
        c := c union {d[i] <= w, -d[i] <= w}
    od;
    ss := simplex[minimize](w, c);
    for i to n + 1 do
        s := op(i, ss);
        if op(1, s) = 'w' then print(s) else assign(s) fi
    od;
    RETURN(evalm(x))
end


  incsolve_weighted := proc(A, b)
local m, n, i, j, c, d, w, x, s, ss;
    if type(A, matrix) = false then
        ERROR(`first parameter must be of type matrix`)
    fi;
    if type(b, vector) = false then
        ERROR(`second parameter must be of type vector`)
    fi;
    m := linalg[rowdim](A);
    n := linalg[coldim](A);
    if m < n then ERROR(`the number of rows must be > or =\
         the number of columns`)
    fi;
    if m <> linalg[vectdim](b) then
        ERROR(`incompatible dimensions`)
    fi;
    x := linalg[vector](n);
    d := linalg[vector](m);
    c := {};
    for i to m do
        d[i] := -b[i];
        for j to n do d[i] := d[i] + A[i, j]*x[j] od;
        c := c union
            {d[i] <= abs(b[i])*w, -d[i] <= abs(b[i])*w}
    od;
    ss := simplex[minimize](w, c);
    for i to n + 1 do
        s := op(i, ss);
        if op(1, s) = 'w' then print(s) else assign(s) fi
    od;
    RETURN(evalm(x))
end

> model_inf := incsolve(A, b);

                           w = .08168983333


               model_inf := [1.023262833, 1.998340333]

> model;

                   y = 2.007287085 x + .9917127333

> yw_coords;

  [Weight(2.953602000, 1), Weight(5.036794000, 1/2),

        Weight(7.004580000, 1/3), Weight(9.098314000, 1/4),

        Weight(11.02952400, 1/5), Weight(12.93161500, 1/6),

        Weight(15.09333500, 1/7), Weight(17.05073800, 1/8),

        Weight(19.07516700, 1/9), Weight(21.04424800, 1/10)]

> model := stats[fit,leastsquare[[x,y],y=a*x + b]]([x_coords,
> yw_coords]);

               model := y = 2.011139185 x + .9705261841

>