read `incsolve.mpl`;
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) + x*rando()/(10.0^5)/10;
od;

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

p1 := plot(points, style=point,color = blue):

A := linalg[matrix](10,2); b := linalg[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];
    yw_coords := [op(y_coords), Weight(y[i],1/i)]; b[i] := y[i];
od;

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

model_inf := incsolve(A, b);

p2 := plot(op(2, model), x=0..10, color = red):

plots[display]({p1, p2});