Two Examples of Least Squares Approximation
<Text-field style="Heading 1" layout="Heading 1">Example 1: Linear Function</Text-field>restart:with(LinearAlgebra):In the manufacture of product Z, the amount of compound A present in the product is controlled by the amount of ingredient B used in the refining process. In manufacturing a litre of Z, the amount of B used and the amount of A present are recorded. The following data were recorded.B:=Matrix(1,5,[[2,4,6,8,10]]);A:=Matrix(1,5,[[3.5,8.2,10.5,12.9,14.6]]);Plot the data. If it seems reasonable after plotting, fit the least squares line of approximation, then replot showing both the line and the points.for i from 1 by 1 to 5 do pts[i]:=[B[1,i],A[1,i]]; end do;with(plots):with(plottools):p1:=plot([seq(pts[i],i=1..5)],x=0..10,y=0..16,style=point,symbol=asterisk,scaling=constrained,color=blue):display(p1);So, it appears that the data is roughly linear. Now we use least squares to compute the best approximation line.M:=Matrix(5,2);ym:=Matrix(5,1);for i from 1 by 1 to 5 do M[i,1]:=1; M[i,2]:=B[1,i]; ym[i,1]:=A[1,i]; end do;M;ym;MtM:=Transpose(M) . M;Mtym:=Transpose(M) . ym;MtMinv:=MatrixInverse(MtM);soln:=MtMinv . Mtym;f:=x->soln[2,1]*x+soln[1,1];f(x);p2:=plot(f(x),x=0..10,y=0..16,scaling=constrained,color=black):display(p1,p2);As we can see, the straight line is pretty good. We can also compute total error.model_data:=M . soln;error_mat:=ym - model_data;? MatrixNormtot_err:=MatrixNorm(error_mat,Euclidean);LUklbXJvd0c2Iy9JK21vZHVsZW5hbWVHNiJJLFR5cGVzZXR0aW5nR0koX3N5c2xpYkdGJzYjLUkjbWlHRiQ2I1EhRic=
<Text-field style="Heading 1" layout="Heading 1">Example 2: Non-Linear Function</Text-field>restart:with(LinearAlgebra):A set of seasonal farm employment data ( employment versus time ) has been collected over a two year period, with time measured in months and employment in millions of people. The belief is that the data should fit a model of the form: f(t) = a + bt + c(cos(t)). Use least squares approximation to determine best fit values for a, b and c.B:=Matrix(1,13,[[3.1,4.3,5.6,6.8,8.1,9.3,10.6,11.8,13.1,14.3,15.6,16.8,18.1]]);A:=Matrix(1,13,[[3.7,4.4,5.3,5.2,4.0,3.6,3.6,5.0,5.0,3.8,2.8,3.3,4.5]]); evalm(A); evalm(B); # B is times, A is employmentPlot the data. If it seems reasonable after plotting, fit the least squares model of approximation, then replot showing both the model and the points.for i from 1 by 1 to 13 do pts[i]:=[B[1,i],A[1,i]]; end do;with(plots):with(plottools):p1:=plot([seq(pts[i],i=1..13)],x=0..24,y=0..6,style=point,symbol=asterisk,scaling=constrained,color=blue):display(p1);So, it appears that the data may follow a cosine type curve. Now we use least squares to compute the best approximation model.M:=Matrix(13,3);ym:=Matrix(13,1);for i from 1 by 1 to 13 do M[i,1]:=1; M[i,2]:=B[1,i]; M[i,3]:=cos(B[1,i]); ym[i,1]:=A[1,i]; end do;evalm(M);evalm(ym);MtM:=Transpose(M) . M;Mtym:=Transpose(M) . ym;MtMinv:=MatrixInverse(MtM);soln:=MtMinv . Mtym;f:=t->soln[1,1]+soln[2,1]*t+soln[3,1]*cos(t);f(t);p2:=plot(f(x),x=0..24,y=0..6,scaling=constrained,color=black):display(p1,p2);As we can see, the model is pretty good. We can also compute total error.model_data:=M . soln;error_mat:=ym - model_data;evalm(error_mat);tot_err:=MatrixNorm(error_mat,Euclidean);JSFH