Linear Fit

points={{1.0,4},{2,0},{3,4},{4,4},{5,6}};
b={4.0,0,4,4,6}; (* b= copy of y-values in points *)
mat={{1.0,1},{2,1},{3,1},{4,1},{5,1}};
(*mat: first column=co[y of x-values in points *)

u=Inverse[Transpose[mat].mat].Transpose[mat].b;
(* linear fit: u[[1]]x+ u[[2]]  *)

plot1=ListPlot[points];
plot2=Plot[  u[[1]]*x + u[[2]],{x,0,5}, PlotStyle-> Red];
Show[plot1,plot2]