Comparison of LMS and LTS Algorithms
/****************************************************************/
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: ROBUST2 */
/* TITLE: Comparison of LMS and LTS Algorithms */
/* PRODUCT: IML */
/* DATA: */
/* */
/* SUPPORT: Rick Wicklin UPDATE: SEP 2013 */
/* REF: */
/* MISC: */
/* */
/****************************************************************/
proc iml;
/* Obs X1 X2 X3 Y Stack Loss data */
SL = { 1 80 27 89 42,
2 80 27 88 37,
3 75 25 90 37,
4 62 24 87 28,
5 62 22 87 18,
6 62 23 87 18,
7 62 24 93 19,
8 62 24 93 20,
9 58 23 87 15,
10 58 18 80 14,
11 58 18 89 14,
12 58 17 88 13,
13 58 18 82 11,
14 58 19 93 12,
15 50 18 89 8,
16 50 18 86 7,
17 50 19 72 8,
18 50 19 79 8,
19 50 20 80 9,
20 56 20 82 15,
21 70 20 91 15 };
x = SL[, 2:4]; y = SL[, 5];
/* Use 2000 Random Subsets for LMS */
optn = j(9,1,.);
optn[2]= 2; /* print a moderate amount of output */
optn[3]= 1; /* compute only LMS regression */
ods select IterHist0 BestSubset EstCoeff;
call lms(sc, coef, wgt, optn, y, x);
ods select all;
r1 = {"Quantile", "Number of Subsets", "Number of Singular Subsets",
"Number of Nonzero Weights", "Min Objective Function",
"Preliminary Scale Estimate", "Final Scale Estimate",
"Robust R Squared", "Asymptotic Consistency Factor"};
sc1 = sc[1:9];
print sc1[r=r1 L="LMS Information and Estimates"];
create Stackloss from SL[c={"ObsNum" "x1" "x2" "x3" "y"}];
append from SL;
close Stackloss;
quit;
proc iml;
use Stackloss;
read all var _num_ into SL;
close Stackloss;
x = SL[, 2:4]; y = SL[, 5];
/* Use 500 random subsets for FAST-LTS algorithm */
optn = j(9,1,.);
optn[2]= 0; /* suppress output */
optn[3]= 0; /* compute only LTS regression */
optn[9]= 0; /* FAST-LTS */
call lts(sc, coef, wgt, optn, y, x);
r1 = {"Quantile", "Number of Subsets", "Number of Singular Subsets",
"Number of Nonzero Weights", "Min Objective Function",
"Preliminary Scale Estimate", "Final Scale Estimate",
"Robust R Squared", "Asymptotic Consistency Factor"};
sc1 = sc[1:9];
print sc1[r=r1 L="LTS Information and Estimates"];
print (coef[1,])[L="Estimated Coefficients"
c={"x1" "x2" "x3" "Intercept"}];
outliers = loc(wgt[1,]=0);
print outliers;