Mathematical Principles for Scientific Computing and Visualization  
http://www.farinhansford.com/books/scv
by Gerald Farin and Dianne Hansford


This notebook illustrates how outliers can create poor results in least squares approximation problem and how we can manage outliers by using weighted least squares approximation.
Figure 8.7 was created with this notebook.

In the first partbelow, we use standard least square approximation to a data set with one outlier. The result is a curve with an unwanted wiggle.

In the second part, we use the weighted technique, which weights each data point based on its distance from the standard least squares solution (previous one).
The result is clearly better. 

See Section 8.4 for details.




See Chapter 8: Data Fitting for details.

Updated July 2008

(* Set the output directory *)
SetDirectory[ToFileName[Extract["FileName"/.NotebookInformation[EvaluationNotebook[]],{1},FrontEnd`FileName]]];



(* Function from which to extract values  at the u_i *)
f[x_]:=N[(Sqrt[x])+((Sin[x])/10)];

(* Function to evaluate a degree n polynomial *)
poly[n_,A_,x_]:=(∑_(j=1)^(n+1))A[[j]] (x^(j-1));

(* Function to solve the normal equations for the polynomial coefficients *)
Coeffs[Vander_,data_]:=Inverse[Transpose[Vander].Vander].(Transpose[Vander].data);


(* Degree polynomial*)
 n=8;

(* Number of data values *)
L=20;

(* Create knots *)
u=Table[i,{i,1,L}];

data=Table[f[u[[i]]],{i,1,L}];

(* Create our outlier *)
 data[[5]]=6;


(* Create the design matrix or Vandermonde matrix for L uniformly space u-values u=1, 2,...L and a degree n polynomial *)
Vander=Table[((u[[i]])^(j-1)),{i,1,L},{j,1,n+1}];


(* Form the matrix part of the normal equations *)
Mat=Transpose[Vander].Vander;

(* Assign the polynomial coefficients to A *)
 A=N[Coeffs[Vander,data]];

(* Create a plot of the data points *)
datplot=ListPlot[data,PlotStyle->{PointSize[0.02`]},BaseStyle->{"FontFamily"->"Helvetica","FontSize"->11}];

(* Create a plot of the given function f and the approximating polynomial *)
PLT=Plot[{f[x],poly[n,A,x]},{x,u[[1]],u[[L]]},PlotStyle->{{Thickness[0.01`]},{Thickness[0.01`]}},BaseStyle->{"FontFamily"->"Helvetica","FontSize"->11}];

(* Put the pics together *)
 pic1=Show[PLT,datplot, PlotRange->All]

 Export["outlier1.eps",pic1] 


(* ========================================= *)
(* Now let's do the weighted least squares *)
(* ========================================= *)


dataold=data; 
(* data[[5]]=6; *)
(* recall that data has outlier data[[5]]=6; *)

(* Compute the distance of the data value to the least squares polynomial from the previous math cell *)

dist=Table[N[poly[n,A,i]-data[[i]]],{i,1,L}];

(* Compute weights for each data *)
(* Uses a heuristic value of 0.1 to prevent zero divide *)
weights=Table[(1/(Abs[dist[[i]]]+0.1`)),{i,1,L}];

Print["weights for each data: ",weights]
Print["  smallest weight should be ",weights[[5]]]

(* The new data are the old ones multiplied by the weights *)
data=Table[weights[[i]] data[[i]],{i,1,L}];

(* Modify the Vandermonde matrix *)
Vander=weights Vander;

(* Recompute the least square approx. *)
A=N[Coeffs[Vander,data]];

(* Create a plot of the data points *)
datplot=ListPlot[dataold,PlotStyle->{PointSize[0.02`]},BaseStyle->{"FontFamily"->"Helvetica","FontSize"->11}];

(* Create a plot of the given function and new approxiamtion *)
PLT=Plot[{f[x],poly[n,A,x]},{x,u[[1]],u[[L]]},PlotStyle->{{Thickness[0.01`]},{Thickness[0.01`]}},BaseStyle->{"FontFamily"->"Helvetica","FontSize"->11}];

 pic2=Show[PLT,datplot, PlotRange->All]

Export["outlier2.eps",pic2]