Chapter Fifteen; Robust Regression:  An Application of Mathematica to Data Analysis

15.1 Fit and Its Problems

lindata = Table[{x, 2 + x + 0.1*Random[ ]}, {x, 0, 5, 0.5}]

quaddata = Table[{x, 2 + x + 0.2 x^2 + 0.1*Random[ ]}, 
{x, 0, 5, 0.5}]

linfit = Fit[lindata, {1, x}, x]

quadfit = Fit[quaddata, {1, x, x^2}, x]

plotfit[fn_, data_] := 
Module[{minind, maxind, mindep, maxdep},
minind = Min[Transpose[data][[1]]]; 
maxind = Max[Transpose[data][[1]]];
mindep = Min[Transpose[data][[2]]]; 
maxdep = Max[Transpose[data][[2]]];
Plot[fn, {x, Floor[minind], Ceiling[maxind]}, 
AxesOrigin -> {0, 0},
PlotRange -> {Min[0, Floor[mindep]], Ceiling[maxdep] + 1},
Evaluate[Epilog -> {PointSize[0.03], Map[Point, data]}]]]

plotfit[linfit, lindata];

plotfit[quadfit, quaddata];

badlindata = lindata;
badlindata[[5]] = {2, 7};
badlindata[[11]] = {5, 1};
badquaddata = quaddata;
badquaddata[[5]] = {2, 15};
badquaddata[[11]] = {5, 1};

badlinfit = Fit[badlindata, {1, x}, x]

badquadfit = Fit[badquaddata, {1, x, x^2}, x]

plotfit[badlinfit, badlindata];

plotfit[badquadfit, badquaddata];

15.2 Least Median of Squares

Median[list_] := Sort[list][[(Length[list]+1)/2]] /;
			OddQ[Length[list]]

Median[list_] := 
	Block[{s, n},
		s = Sort[list] ;
		n = Length[list] ;
		(s[[n/2]] + s[[n/2 + 1]]) / 2
	] /; EvenQ[Length[list]]

15.2.1 Conceptual Approach: 2 Dimensions

lmsfita[data_, x_] := Module[
{l, v, i, j, w, h, a, sl, dx},
l = Length[data];
v = Flatten[Table[{data[[i]], data[[j]]}, {i,1,l-1}, {j,i+1,l}], 1];
w = Union[Map[
 Module[{dx},
 dx = #[[1, 1]] - #[[2, 1]]; 
 If[dx !=0, 
 {(#[[1, 2]]-#[[2, 2]])/dx, 
 (#[[1, 1]] #[[2, 2]] - #[[1, 2]] #[[2, 1]])/dx},
 {Infinity, #[[1, 1]]}]
 ]&, 
	v]];
h = Function[a, 
{Median[Map[ 
(If[a[[1]] != Infinity, 
(#[[2]] - a[[1]] #[[1]] - a[[2]])^2,
(a[[2]] - #[[1]])^2])&, data]],
 a[[1]], a[[2]]}];
sl = Sort[Map[h, w]][[1]];
If[sl[[2]] != Infinity, sl[[2]] x + sl[[3]], 
Print["x = ",sl[[3]]]]
]

DefineTheLines[data_] := Module[
{l, v, i, j},
l = Length[data];
v = Flatten[Table[{data[[i]], data[[j]]}, 
		{i,1,l-1}, {j,i+1,l}], 1];
v]

DefineTheLines[{{1, 2}, {2, 3}, {4, 5}}]

lmslinfit = lmsfita[badlindata, x]

plotfit[lmslinfit, badlindata];

vertdata = {{1, 1}, {1, 2}, {1, 3}, {4, 5}};

lmslinfit = lmsfita[vertdata, x]

15.3 Higher Dimensions

KSubsets[l_List,0] := { {} }
KSubsets[l_List,1] := Partition[l,1]
KSubsets[l_List,k_Integer?Positive] := {l} /; (k == Length[l])
KSubsets[l_List,k_Integer?Positive] := {}  /; (k > Length[l])
KSubsets[l_List,k_Integer?Positive] :=
	Join[
		Map[(Prepend[#,First[l]])&, KSubsets[Rest[l],k-1]],
		KSubsets[Rest[l],k]
		]

KSubsets[{1, 2, 4, 3}, 1]

Vertices = Table[{Cos[2 Pi i/20], Sin[2 Pi i/20]}, {i, 1, 20}];
Show[Graphics[{
{Thickness[0.002], RGBColor[0, 0, 1], Map[Line, KSubsets[Vertices, 2]]}, 
{PointSize[0.03], RGBColor[1, 0, 0], Map[Point, Vertices]}}, 
AspectRatio -> 1]];

RandomPermutation[n_Integer?Positive] :=
Map[Last, Sort[Map[({Random[], #})&, Range[n]]]]

myperm = RandomPermutation[5]

mylist = {a, b, c, d, e};
mylist[[myperm]]

Unprotect[Table];
Table[expr_, {l__List}] := Table[expr, l];
Protect[Table];

PSubsets[l_List, n_Integer?Positive] :=
Module[{itlist, expr,a},
a[0] = 0;
itlist = 
     Table[{a[j], a[j-1] + 1, Length[l] + j - n}, {j, 1, n}];
expr = Table[l[[Array[a, {n}]]], Evaluate[itlist]];
expr = Flatten[expr, n-1]
] /; ( n <= Length[l])

mylist = {a, b, c, d, e};
PSubsets[mylist, 2]

averes[data_] := Sum[Timing[KSubsets[data, n];][[1]], {n, 2, 5}]/
                 Sum[Timing[PSubsets[data, n];][[1]], {n, 2, 5}];

{averes[Range[10]], averes[Range[20]]}

Fitlms[data_List, fns_List, vbls_List, thoroughness_:1] :=
Module[{makelinear, lineardata, vlen, flen, calcmedian,
aa, hyperplane, hyperplanes, numplanes, eqn, solver, solns, 
coeffs, sl},
vlen = Length[vbls]; flen = Length[fns];

(* makelinear takes a function list like {1,x,x y}, together with a
   variable list --- in this instance {x,y} --- to create the function
   {1,#1, #1 #2, #3}& *)

makelinear = Evaluate[Join[
fns /. Table[vbls[[i]] -> Slot[i], {i, 1, vlen}],
  {Slot[vlen + 1]}]
		]&;

(* now we generate all hyperplane, see how many there are and take a 
   random subset if thoroughness < 1 *)

lineardata = Apply[makelinear, data, 1];
hyperplanes = KSubsets[lineardata, flen];
numplanes = Length[hyperplanes];
hyperplanes = If[thoroughness < 1, 
Take[hyperplanes[[RandomPermutation[numplanes]]], 
Ceiling[thoroughness*numplanes]],
 hyperplanes];

(* eqn takes a hyperplane and constructs the relevant linear equations to
   determine the associated values of the parameters  *)

coeffs = Array[aa, flen];
eqn = Function[hyperplane,
Dot[coeffs, Drop[hyperplane, -1]] == 
Last[hyperplane]];

(* Now construct the parameters for all the hyperplanes chosen above *)

solver = coeffs /. Solve[Map[eqn, #], coeffs]&;
solns = Flatten[Map[solver, hyperplanes], 1];

(* finally for each hyperplane, the median-square-deviation is computed, 
   and the hyperplane for which this is a minimum is returned *)

calcmedian = 
Function[coeffset, 
{Median[Map[(Last[#] - Dot[coeffset, Drop[#, -1]])^2&,lineardata]], 
coeffset}];
sl = Sort[Map[calcmedian, solns]][[1]];
Dot[sl[[2]], fns]
]

lms = Fitlms[badquaddata, {1, x, x^2}, {x}, 0.2]

plotfit[lms, badquaddata];

RandomKSubset[n_Integer,k_Integer] := RandomKSubset[Range[n],k]

RandomKSubset[set_List,k_Integer] := 
 Module[{s=Range[Length[set]],i,n=Length[set],x},
 set [[
 Sort[
 Table[
 x=Random[Integer,{1,i}];
 {s[[i]],s[[x]]} = {s[[x]],s[[i]]};
 s[[i]],
 {i,n,n-k+1,-1}
 ]
 ]
 ]]
]

Fitlms[data_List, fns_List, vbls_List, 
       thoroughness_:1, method_Integer:1] :=
Module[{makelinear, lineardata, vlen, flen, calcmedian,
aa, hyperplane, hyperplanes, numplanes, eqn, solver, solns, 
coeffs, sl},
vlen = Length[vbls]; flen = Length[fns];
makelinear = 
Evaluate[Join[ fns /. Table[vbls[[i]] -> Slot[i], {i, 1, vlen}],
		{Slot[vlen + 1]}]
		]&;
lineardata = Apply[makelinear, data, 1];
datano = Length[lineardata];
numplanes = datano!/(flen! (datano-flen)!);
If[method == 1,
     hyperplanes = KSubsets[lineardata, flen];
     hyperplanes = If[thoroughness < 1, 
     Take[hyperplanes[[RandomPermutation[numplanes]]], 
     Ceiling[thoroughness*numplanes]],
             hyperplanes],

     (* method 2 *)
     hyperplanes = Table[RandomKSubset[lineardata,flen],
                               {thoroughness*numplanes}]
   ];
coeffs = Array[aa, flen];
eqn = Function[hyperplane,
Dot[coeffs, Drop[hyperplane, -1]] == 
Last[hyperplane]];
solver = coeffs /. Solve[Map[eqn, #], coeffs]&;
solns = Flatten[Map[solver, hyperplanes], 1];
calcmedian = 
Function[coeffset, 
{Median[Map[(Last[#] - Dot[coeffset, Drop[#, -1]])^2&,lineardata]], 
coeffset}];
sl = Sort[Map[calcmedian, solns]][[1]];
Dot[sl[[2]], fns]
]

Timing[Fitlms[badquaddata, {1, x, x^2}, {x}, 0.2,2]]

Timing[Fitlms[badquaddata,{1,x,x^2},{x},0.2,1]]

Nothing in it! What about a larger data set?

big = Table[{x, 2 + x + 0.2 x^2 + 0.1*Random[ ]}, 
{x, 0, 15, 0.5}];

Timing[Fitlms[big, {1, x, x^2}, {x}, 0.01,2]]

Timing[Fitlms[big, {1, x, x^2}, {x}, 0.01,1]]