(* Robust regression methods *)

(* basic stuff *)
MEDIAN[x_]:=Module[{X=x,L,M},
L=Length[X];X=Sort[X];
If[Mod[L,2]==0,M=0.5(X[[L/2]]+X[[L/2+1]]),M=X[[(L+1)/2]]];
M];
MAD[x_]:=Module[{X=x,A},
A=Abs[X-MEDIAN[X]];
MEDIAN[A]];
MADN[x_]:=MAD[x]/0.6745;

(* bisquare functions for iterative reweighted regression *)
δ=1/2;
WBisquare[t_]:=Min[3-3(t^2)+(t^4),1/(t^2)];
wBisquare[x_,σ_]:=WBisquare[x/σ];
sigBisquare[x_,σ_]:=Sqrt[Sum[wBisquare[x[[i]],σ]((x[[i]])^2),{i,1,Length[x]}]/Length[x]/δ];

(* data *)
a=1;
b=2;
(* first part of data *)
x=RandomReal[UniformDistribution[{1,4}],50];
e=RandomReal[NormalDistribution[0,0.2],50];
y=a x+b+e;
(* contamination with outliers *)
xcont=RandomReal[NormalDistribution[7,0.5],30];
ycont=RandomReal[NormalDistribution[2,0.5],30];
x=Join[x,xcont];
y=Join[y,ycont];
n=Length[x];
w=Table[0,{i,1,n}];
data=Table[{x[[i]],y[[i]]},{i,1,n}];

(* standard least-squares regression *)
line=Fit[data,{1,X},X];
a0LS=SeriesCoefficient[Series[line,{X,0,0}],0];
a1LS=SeriesCoefficient[Series[line,{X,1,1}],1];
residuals0=Table[y[[i]]-(a0LS+a1LS x[[i]]),{i,1,n}];
sighat=MEDIAN[Abs[residuals0]]/0.6745;
Print["standard regression: ",a0LS," ",a1LS," ",sighat];

(* iterative reweighted LS regression *)
(* start with regular LS *)
a0IRWLS=a0LS;
a1IRWLS=a1LS;
For[k=1,k<30,k++,{
(* residuals *)
r=Table[y[[i]]-(a0IRWLS+a1IRWLS x[[i]]),{i,1,n}];
For[i=1,i≤n,i++,w[[i]]=wBisquare[r[[i]],sighat]];
A=Sum[w[[i]],{i,1,n}];
Ap=B=Sum[w[[i]]x[[i]],{i,1,n}];
Bp=Sum[w[[i]]((x[[i]])^2),{i,1,n}];
CC=Sum[w[[i]]y[[i]],{i,1,n}];
CCp=Sum[w[[i]]x[[i]]y[[i]],{i,1,n}];
sol=Solve[{aa0 A+aa1 B==CC,aa0 Ap+aa1 Bp==CCp},{aa0,aa1}];
{{a0IRWLS,a1IRWLS}}={aa0,aa1}/.sol;
sighat=Sqrt[Sum[w[[i]]((r[[i]])^2),{i,1n}]/n/δ];
}];
Print["iterative reweighted: ",a0IRWLS," ",a1IRWLS," ",sighat];

(* LMS: Least Median Squares *)
a0LMS=0;
LMSFUNC[a1_?NumericQ]:=Module[{A1=a1},
sol=FindMinimum[MEDIAN[(((y- A1 x-A0))^2)],{A0,2}];
a0LMS=A0/.sol[[2]];
MEDIAN[(((y- A1 x-a0LMS))^2)]
];
Plot[LMSFUNC[A1],{A1,0.8,1.2}]
sol=NMinimize[LMSFUNC[A1],{A1}] ;
a1LMS=A1/.sol[[2]];
Print["Least Median Squares: ",a0LMS," ",a1LMS];
Plot3D[Log[MEDIAN[(((y-a1LMS x-a0LMS))^2)]],{a0LMS,1.7,2.2},{a1LMS,0.7,1.3},PlotPoints->50]
ContourPlot[Log[MEDIAN[(((y-a1LMS x-a0LMS))^2)]],{a0LMS,1.7,2.2},{a1LMS,0.7,1.3},Contours->50,ColorFunction->"Pastel"]

Show[ListPlot[data,PlotRange->{{0,10},{0,15}}],Plot[line,{X,0,10},PlotStyle->Blue],Plot[a0IRWLS+a1IRWLS X,{X,0,10},PlotStyle->Red],Plot[a0LMS+a1LMS X,{X,0,10},PlotStyle->Green]]

(* output to files *)

out=OpenWrite["ObjectiveFunc.dat"];
For[AA0=1.7,AA0≤2.2001,AA0+=0.017241379,{
For[AA1=0.85,AA1≤1.15001,AA1+=0.010344828,{
OF=Log[MEDIAN[(((y-AA1 x-AA0))^2)]];
WriteString[out,ToString[CForm[SetPrecision[OF,7]]]," "];
}]
WriteString[out,"\n"];
}]
Close[out];

out=OpenWrite["LMSdata.dat"];
For[i=1,i≤n,i++,{
WriteString[out,ToString[CForm[SetPrecision[x[[i]],7]]]," ",ToString[CForm[SetPrecision[y[[i]],7]]],"\n"];
}]
Close[out];