Fitting a polynomial of degree n-1 to n points

p1 = {0,3};
p2 = {1,1};
p3 = { 3,4};
p4 = {6,2};
p = {p1,p2,p3,p4};

points = 
ListPlot[
p,
PlotStyle->Directive[Large,Red],
AxesStyle->Thick,
TicksStyle->Large,
PlotRange->{{-1,7},All},
ImageSize->720
]

Vmat[p_] :=
Table[
Which[
i==1&&j==0,1,
True,p[[i,1]]^j
],
{i,1,Length[p]},{j,0,Length[p]-1}
];
Vmat[p]//MatrixForm

bvec[p_] := 
Table[p[[i,2]],{i,1,Length[p]}];
bvec[p]

coefs = LinearSolve[Vmat[p],bvec[p]]

poly[x_] := Sum[coefs[[i]]*x^(i-1),{i,1,Length[coefs]}];

polyplot =
Plot[
poly[x],{x,-1,7},
PlotStyle->Directive[Blue,Thick],
AxesStyle->Thick,
TicksStyle->Large,
PlotRange->All,
ImageSize->720
];
Show[polyplot,points]

Fitting a polynomial of degree n-1 to n points by hand

Newtoncoefs = {3,-2,7/6,-4/15};
Newtonpoly[x_] :=
Sum[Newtoncoefs[[i]]*Product[x-p[[j,1]],{j,1,i-1}],{i,1,Length[Newtoncoefs]}];

Newtonpoly[x]

poly[x]

Expand[Simplify[Newtonpoly[x]]]

Lagrange’s formula for a polynomial of degree n-1 through n points

lagrangepoly[i_,x_] :=
Product[
Which[
j≠i,x-p[[j,1]],
True,1
],{j,1,Length[p]}]/
Product[
Which[
j≠i,p[[i,1]]-p[[j,1]],
True,1
],{j,1,Length[p]}];

l1 = 
Plot[
p[[1,2]]*lagrangepoly[1,x],{x,-1,7},
PlotStyle->Directive[Thick,Orange],
AxesStyle->Thick,
TicksStyle->Large,
PlotRange->All,
ImageSize->720
];
l2 =
Plot[
p[[2,2]]*lagrangepoly[2,x],{x,-1,7},
PlotStyle->Directive[Thick,Yellow],
AxesStyle->Thick,
TicksStyle->Large,
PlotRange->All,
ImageSize->720
];
l3 =
Plot[
p[[3,2]]*lagrangepoly[3,x],{x,-1,7},
PlotStyle->Directive[Thick,Darker[Green]],
AxesStyle->Thick,
TicksStyle->Large,
PlotRange->All,
ImageSize->720
];
l4 =
Plot[
p[[4,2]]*lagrangepoly[4,x],{x,-1,7},
PlotStyle->Directive[Thick,Purple],
AxesStyle->Thick,
TicksStyle->Large,
PlotRange->All,
ImageSize->720
];

Show[
l1,l2,l3,l4,points
]

Overfitting example

m = 1;
b = 1;

p1 = {0,m*0+b + RandomVariate[NormalDistribution[0,0.5]]};
p2 = {1,m*1+b + RandomVariate[NormalDistribution[0,0.5]]};
p3 = {3,m*3+b + RandomVariate[NormalDistribution[0,0.5]]};
p4 = {6,m*6+b + RandomVariate[NormalDistribution[0,0.5]]};
p = {p1,p2,p3,p4};

line =
Plot[
m*x + b,{x,-1,7},
PlotStyle->Directive[Large,Darker[Green]],
AxesStyle->Thick,
TicksStyle->Large,
PlotRange->{{-1,7},All},
ImageSize->720
];
points = 
ListPlot[
p,
PlotStyle->Directive[Large,Red],
AxesStyle->Thick,
TicksStyle->Large,
PlotRange->{{-1,7},All},
ImageSize->720
];

Show[
line,points
]

Vmat[p_] :=
Table[
Which[
i==1&&j==0,1,
True,p[[i,1]]^j
],
{i,1,Length[p]},{j,0,Length[p]-1}
];
Vmat[p]//MatrixForm

bvec[p_] := 
Table[p[[i,2]],{i,1,Length[p]}];
bvec[p]

coefs = LinearSolve[Vmat[p],bvec[p]]
poly[x_] := Sum[coefs[[i]]*x^(i-1),{i,1,Length[coefs]}];

polyplot =
Plot[
poly[x],{x,-1,7},
PlotStyle->Directive[Blue,Thick],
AxesStyle->Thick,
TicksStyle->Large,
PlotRange->All,
ImageSize->720
];
Show[polyplot,line,points]