SAS/IML V8 Modules for Robust Regression
/****************************************************************/
/* S A S S A M P L E L I B R A R Y */
/* */
/* NAME: ROBUSTMC */
/* TITLE: SAS/IML V8 Modules for Robust Regression */
/* PRODUCT: IML */
/* SYSTEM: ALL */
/* KEYS: */
/* PROCS: */
/* DATA: */
/* */
/* SUPPORT: GJW UPDATE: FEB 97 */
/* REF: */
/* MISC: */
/* Modules: PRILMTS - PRINTING LMS AND LMT RESULTS
/* SCATLMS - PLOTTING LMS AND LMT RESULTS
/* PRIMVE - PRINTING MVE RESULTS
/* SCATMVE - PLOTTING MVE RESULTS
/* SCATMCD - PLOTTING MCD RESULTS
/* LMSDIAP - DIAGNOSTIC OUTLIER PLOT
/****************************************************************/
proc iml;
/*------------------------------------------------------------*/
start prilmts(job,sc,coef,wgt);
/* job=1: LMS
=2: LMS + LS
=3: LMS + RLS
=4: LTS
=5: LTS + LS
=6: LTS + RLS */
if job < 4 then do;
title1= "Results of Least Median Squares Estimation";
title2= "Estimated LMS Coefficients";
title3= "LMS Residuals";
end; else do;
title1= "Results of Least Trimmed Squares Estimation";
title2= "Estimated LTS Coefficients";
title3= "LTS Residuals";
end;
print title1[label=""];
print "Quantile. . . . . . . . . . ." (sc[1]),
"Number of Subsets. . . . . . " (sc[2]),
"Number of Singular Subsets . " (sc[3]),
"Number of Nonzero Weights. . " (sc[4]),
"Objective Function. . . . . ." (sc[5]),
"Preliminary Scale Estimate. ." (sc[6]),
"Final Scale Estimate. . . . ." (sc[7]),
"Robust R Squared. . . . . . ." (sc[8]),
"Asymptotic Consistency Factor" (sc[9]);
if job = 2 | job = 5 then do;
print "LS Scale Estimate . . . . ." (sc[11]),
"Sum of Squares. . . . . . ." (sc[12]),
"R-squared . . . . . . . . ." (sc[13]),
"F Statistic . . . . . . . ." (sc[14]);
end;
if job = 3 | job = 6 then do;
print "RLS Scale Estimate. . . . ." (sc[11]),
"Weighted Sum of Squares . ." (sc[12]),
"Weighted R-squared. . . . ." (sc[13]),
"F Statistic . . . . . . . ." (sc[14]);
end;
nr2 = nrow(coef);
if nr2 > 0 then do;
print (coef[1,])[label=title2];
print (coef[2,])[label="Indices of Best Sample"];
if nr2 > 2 then do;
print (coef[3,])[label="Estimated WLS Coefficients"];
if nr2 = 8 then do;
print (coef[4,])[label="Standard Errors"];
print (coef[5,])[label="T Values"];
print (coef[6,])[label="Probabilities"];
print (coef[7,])[label="Lower Wald CI"];
print (coef[8,])[label="Upper Wald CI"];
end; end; end;
nr3 = nrow(wgt);
if nr3 > 0 then do;
print (wgt[1,])[label="Estimated Weights"];
print (wgt[2,])[label=title3];
if nr3 > 2 then do;
print (wgt[3,])[label="Diagnostics"];
if nr3 = 4 then do;
print (wgt[4,])[label="WLS Residuals"];
end; end; end;
finish prilmts;
/*------------------------------------------------------------*/
start scatlms(job,optn,x,y,xnam,ynam,titl,ipsf,filn);
/* job = 0: only draw scatter plot without regression
= 1: LS
= 2: LS + LMS
= 3: LS + LMS + WLS
= 4: LS + LTS
= 5: LS + LTS + WLS
= 6: LS + LMS + LTS
= 7: LS + LMS + LTS + WLS (WLS same LMS and LTS)
* the next 8-14 are 1-7 but with L1 regression *
= 8: LS + L1
= 9: LS + L1 + LMS
=10: LS + L1 + LMS + WLS
=11: LS + L1 + LTS
=12: LS + L1 + LTS + WLS
=13: LS + L1 + LMS + LTS
=14: LS + L1 + LMS + LTS + WLS (WLS same LMS and LTS)
ipsf: specifies GOPTIONS statement:
ipsf=0: draw interactively
ipsf=1: write file with smaller picture for document
ipsf=2: write file with larger picture for color printer
filn: name for graphic segment */
start gscenter(x,y,str);
call gstrlen(len,str);
call gscript(x-len/2,y,str);
finish gscenter;
start goptlms(ipsf,filn);
if ipsf then do;
fnam = (filn || ".ps");
fnam = rowcatc(fnam);
print fnam;
call execute("filename gsasfile'",fnam,"';");
if ipsf = 1 then do;
call execute('goptions reset=goptions dev=pslepsf
gsfmode=replace gaccess=gsasfile gsfname=GSASFILE
hsize=5.625 in vsize=4.65 in htext=3.0 pct border;');
end; else do;
call execute('goptions reset=goptions dev=pscolor
gsfmode=replace gaccess=gsasfile border erase;');
end;
end; else do;
call execute('goptions reset=goptions /*dev=xcolor*/;');
end;
finish goptlms;
noint = optn[1]; nobs = nrow(x); n = ncol(x);
if nrow(y) ^= nobs | ncol(y) ^= 1 | n ^= 1 then do;
print "Data vectors x and y not compatible";
goto leave;
end;
/* get [xmin,xmax] and [ymin,ymax] */
xmin1 = x[><]; xmax1 = x[<>];
xd1 = xmax1 - xmin1; xc1 = .1 * xd1;
xmin2 = xmin1 - xc1; xmax2 = xmax1 + xc1;
xd2 = xmax2 - xmin2; xc2 = .1 * xd2;
xme = .5*(xmin2+xmax2);
ymin1 = y[><]; ymax1 = y[<>];
yd1 = ymax1 - ymin1;
ils = 0; iwls = 0; ilms = 0; ilts = 0; il1 = 0;
if job > 0 then ils = 1;
if job = 3 | job = 5 | job = 7 then iwls = 1;
if job =10 | job =12 | job =14 then iwls = 1;
if job = 2 | job = 3 | job = 6 | job = 7 then ilms = 1;
if job = 9 | job =10 | job =13 | job =14 then ilms = 1;
if job = 4 | job = 5 | job = 6 | job = 7 then ilts = 1;
if job =11 | job =12 | job =13 | job =14 then ilts = 1;
if job > 7 & job < 15 then il1 = 1;
/* print "LS=" ils "WLS=" iwls "L1=" il1; */
/* print "LMS=" ilms "LTS=" ilts; */
/* enumerate outlier points if either LMS or LTS */
klmts = 0;
if job = 2 | job = 3 | job = 9 | job =10 then klmts = 1;
if job = 4 | job = 5 | job =11 | job =12 then klmts = 2;
nr = 0;
kls = 0; kwls = 0; klms = 0; klts = 0; kl1 = 0;
vreg = j(10,2,0); nreg = j(10,1," ");
if job then do;
/* A. Do LMS and LTS */
if (ilms) then do;
/*--- 1. Compute LS and LMS ---*/
optn[3]= 2; /* ilsq: LS and LMS */
/* print "LS and LMS"; */
CALL LMS(sc,coef,wgt,optn,y,x);
/* print "End of LMS" coef; */
kls= 1; b = coef[3,1];
if noint=0 then a = coef[3,2]; else a = 0.;
nr = nr + 1; nreg[nr]= "LS";
vreg[nr,1]= xmin1 * b + a;
vreg[nr,2]= xmax1 * b + a;
klms=1; b = coef[1,1];
if noint=0 then a = coef[1,2]; else a = 0.;
nr = nr + 1; nreg[nr]= "LMS";
vreg[nr,1]= xmin1 * b + a;
vreg[nr,2]= xmax1 * b + a;
if ilts then do;
/*--- 2. Compute LTS and WLS ---*/
if iwls then optn[3]= 1; /* ilsq: LTS and WLS */
else optn[3]= 0;
/* print "LTS and WLS"; */
CALL LTS(sc,coef,wgt,optn,y,x);
/* print "End of LTS" coef; */
klts=1; b = coef[1,1];
if noint=0 then a = coef[1,2]; else a = 0.;
nr = nr + 1; nreg[nr]= "LTS";
vreg[nr,1]= xmin1 * b + a;
vreg[nr,2]= xmax1 * b + a;
if iwls then do;
kwls=1; b = coef[3,1];
if noint=0 then a = coef[3,2]; else a = 0.;
nr = nr + 1; nreg[nr]= "WLS";
vreg[nr,1]= xmin1 * b + a;
vreg[nr,2]= xmax1 * b + a;
end;
end; else if iwls then do;
/*--- 3. Compute LMS and WLS: if no LTS is asked ---*/
optn[3]= 1; /* ilsq: LMS and WLS */
/* print "LMS"; */
CALL LMS(sc,coef,wgt,optn,y,x);
/* print "End of LMS" coef; */
kwls=1; b = coef[3,1];
if noint=0 then a = coef[3,2]; else a = 0.;
nr = nr + 1; nreg[nr]= "WLS";
vreg[nr,1]= xmin1 * b + a;
vreg[nr,2]= xmax1 * b + a;
end;
end; else if ilts then do;
/* 4. Do LS and LTS: if no LMS is asked */
optn[3]= 2; /* ilsq: LS and LTS */
/* print "LS and LTS"; */
CALL LTS(sc,coef,wgt,optn,y,x);
/* print "End of LTS" coef; */
kls= 1; b = coef[3,1];
if noint=0 then a = coef[3,2]; else a = 0.;
nr = nr + 1; nreg[nr]= "LS";
vreg[nr,1]= xmin1 * b + a;
vreg[nr,2]= xmax1 * b + a;
klts=1; b = coef[1,1];
if noint=0 then a = coef[1,2]; else a = 0.;
nr = nr + 1; nreg[nr]= "LTS";
vreg[nr,1]= xmin1 * b + a;
vreg[nr,2]= xmax1 * b + a;
if iwls then do;
/*--- 3. Compute LTS and WLS: if no LMS is asked ---*/
optn[3]= 1; /* ilsq: LMS and WLS */
CALL LTS(sc,coef,wgt,optn,y,x);
kwls=1; b = coef[3,1];
if noint=0 then a = coef[3,2]; else a = 0.;
nr = nr + 1; nreg[nr]= "WLS";
vreg[nr,1]= xmin1 * b + a;
vreg[nr,2]= xmax1 * b + a;
end;
end;
/* B. Do L1 and leftover LS */
if il1 then do;
l1opt = j(2,1,0); l1opt[1]=1000;
if noint=0 then xx = x || j(nobs,1,1.);
else xx = x;
CALL LAV(rc,coef,xx,y,,l1opt);
kl1=1; b = coef[1,1];
if noint=0 then a = coef[1,2]; else a = 0.;
nr = nr + 1; nreg[nr]= "L1";
vreg[nr,1]= xmin1 * b + a;
vreg[nr,2]= xmax1 * b + a;
end;
if ils && kls=0 then do;
if noint=0 then do;
xx = x || j(nobs,1,1.);
sxx = xx` * xx;
sxy = xx` * y;
xls = inv(sxx) * sxy;
kls=1; b = xls[1]; a = xls[2];
end; else do;
sxx = x` * x; sxy = x` * y;
kls=1; b = sxy / sxx; a = 0.;
end;
nr = nr + 1; nreg[nr]= "LS";
vreg[nr,1]= xmin1 * b + a;
vreg[nr,2]= xmax1 * b + a;
end;
/* correct ymin and ymax */
do ir=1 to nr;
a = vreg[ir,1]; b = vreg[ir,2];
if a < ymin1 then ymin1 = a;
if b < ymin1 then ymin1 = b;
if a > ymax1 then ymax1 = a;
if b > ymax1 then ymax1 = b;
end;
end;
/* print "NR=" nr; */
/* now with corrected ymin and ymax */
yd1 = ymax1 - ymin1; yc1 = .1 * yd1;
ymin2 = ymin1 - yc1; ymax2 = ymax1 + yc1;
yd2 = ymax2 - ymin2; yc2 = .1 * yd2;
yme = .5*(ymin2+ymax2);
/* No Intercept: regression through origin;
start axis with zero */
if noint^=0 & xmin1 > 0 & xmin2 < xd2
& ymin2 < yd2 then do;
xmin1 = 0.; xmin2 = 0.;
xd2 = xmax2 - xmin2; xc2 = .1 * xd2;
xme = .5*(xmin2+xmax2);
ymin1 = 0.; ymin2 = 0.;
yd2 = ymax2 - ymin2; yc2 = .1 * yd2;
yme = .5*(ymin2+ymax2);
do ir=1 to nr;
vreg[ir,1] = 0;
end;
end;
/* prepare the plot */
xbox = { 0 100 100 0};
ybox = { 0 0 100 100};
wind = (xmin2 || ymin2) // (xmax2 || ymax2);
/* create the plot */
call goptlms(ipsf,filn);
call gstart;
call gopen(filn);
call gset("font","swiss"); /* set character font */
call gpoly(xbox,ybox); /* draw a box around plot */
call gwindow(wind); /* define window */
call gport({15 15 , 85 85}); /* define viewport */
call gxaxis((xmin2 || ymin2),xd2,11, , ,"4.",1.);
call gyaxis((xmin2 || ymin2),yd2,11, , ,"4.",1.);
/* plot xnam, ynam variable names of axis and the title */
call gset("height",1.5);
call gstrlen(len,xnam);
call gtext(xmax2-len,ymin2-yc2,xnam);
call gvtext(xmin2-xc2,ymax2,ynam);
call gset("height",2.0);
call gscenter(xme,ymax2,titl);
/* draw regression lines */
xa = xmin1; xb = xmax1;
do ir=1 to nr;
ya = vreg[ir,1]; yb = vreg[ir,2];
call gdraw((xa || xb),(ya || yb));
call gscript(xb,yb,nreg[ir]);
end;
/* print the scatter points */
call gpoint(x,y,"dot","magenta",1.5);
/* enumerate the zero weight points */
if klmts>0 then do;
do i=1 to nobs;
if wgt[i] = 0 then do;
z = char(i,2,0);
nzwx = x[i] + .01 * xd2;
nzwy = y[i];
call gscript(nzwx,nzwy,z,,,1.);
end; end; end;
call gshow;
call gclose;
call gstop;
leave:
finish scatlms;
/*------------------------------------------------------------*/
start primve(sc,xmve,dist);
n = ncol(xmve);
nr = nrow(xmve);
print "Quantile . . . . . . . . ." (sc[1]),
"Number of Subsets . . . . " (sc[2]),
"Number of Singular Subsets" (sc[3]),
"Number of Nonzero Weights " (sc[4]),
"Minimum Volume . . . . . ." (sc[5]),
"Distance Quantile. . . . ." (sc[6]),
"Cutoff Value . . . . . . ." (sc[7]),
"Correction Factor. . . . ." (sc[8]),
"Scaling Factor . . . . . ." (sc[9]);
print "Robust Location", (xmve[1,]);
print "Robust Scatter Matrix", (xmve[3:2+n,]);
print "Eigenvalues of Robust Scatter Matrix", (xmve[2,]);
if (nr > n+2) then
print "Robust Correlation Matrix", (xmve[3+n:2+n+n,]);
print "Mahalanobis Distances", (dist[1,]);
print "Robust Distances", (dist[2,]);
print "Estimated Weights", (dist[3,]);
finish primve;
/*------------------------------------------------------------*/
start scatmve(job,optn,prob,a,vnam,titl,ipsf,filn);
/* job=1: draw confidence ellipse for mean,
job=2: draw confidence ellipse for prediction
ipsf: specifies GOPTIONS statement:
ipsf=0: draw interactively
ipsf=1: write file with smaller picture for document
ipsf=2: write file with larger picture for color printer
filn: prefix name of postscript file */
start gscenter(x,y,str);
call gstrlen(len,str);
call gscript(x-len/2,y,str);
finish gscenter;
start goptmve(ipsf,filn,n,jx,jy);
if ipsf then do;
if n < 10 then do;
post = jy * 10 + jx;
post = char(post,2,0);
end; else do;
post = jy * 100 + jx;
post = char(post,4,0);
end;
fnam = (filn || post || ".ps");
fnam = rowcatc(fnam);
print fnam;
call execute("filename gsasfile'",fnam,"';");
if ipsf = 1 then do;
call execute('goptions reset=goptions dev=pslepsf
gsfmode=replace gaccess=gsasfile gsfname=GSASFILE
hsize=5.625 in vsize=4.65 in htext=3.0 pct border;');
end; else do;
call execute('goptions reset=goptions dev=pscolor
gsfmode=replace gaccess=gsasfile border erase;');
end;
end; else do;
call execute('goptions reset=goptions /*dev=xcolor*/;');
end;
finish goptmve;
start elli1(job,prob,nobs,xm,covm,nn,corn,elp);
/* This is my version */
/* job=1: draw confidence ellipse for mean,
job=2: draw confidence ellipse for prediction */
tf = 2. * (nobs-1) / (nobs * (nobs - 2));
ff = tf * FINV(prob,2,nobs-2);
if job = 2 then ff = ff * (nobs + 1.);
ff = sqrt(ff);
/* Eigenvectors are stored in columns */
/* get orthogonal rotation matrix */
/* since length of eigenvectors is 1:
rot[1,1]= rot[2,2]= cos(alfa),
rot[2,1]= sin(alfa) = -rot[1,2] */
CALL EIGEN(eval,evec,covm);
if covm[1,1] > covm[2,2] then do;
t = eval[1]; eval[1] = eval[2]; eval[2] = t;
t = evec[1,1]; evec[1,1] = evec[1,2]; evec[1,2] = t;
t = evec[2,1]; evec[2,1] = evec[2,2]; evec[2,2] = t;
end;
if eval[1] <= 0 then da = 0.;
else da = sqrt(eval[1]);
if eval[2] <= 0 then db = 0.;
else db = sqrt(eval[2]);
rot = j(2,2,0.);
rot[1,1] = da * ff * evec[1,1];
rot[2,1] = da * ff * evec[1,2];
rot[1,2] = -db * ff * evec[1,2];
rot[2,2] = db * ff * evec[1,1];
/* Get axes of ellipsoid: a axis: 1 -> 2
b axis: 3 -> 4 */
corn = j(4,2,0.);
corn[1,1] = -1.; corn[1,2] = 0.;
corn[2,1] = 1.; corn[2,2] = 0.;
corn[3,1] = 0.; corn[3,2] = -1.;
corn[4,1] = 0.; corn[4,2] = 1.;
corn = corn * rot` + j(4,1,1) * xm;
/* Points of ellipsoid: elp[nn+1,2] */
pi = 3.141592653589; pi2 = 2. * pi;
elp = j(nn+1,2,0.);
phi = 0.; dp = pi2 / nn;
do j = 1 to nn+1;
x = cos(phi); y = sin(phi);
elp[j,] = (x || y) * rot` + xm;
phi = phi + dp;
end;
finish elli1;
start elli2(job,prob,nobs,xm,covm,nn,corn,elp);
/* This is the INSIGHT version */
/* job=1: draw confidence ellipse for mean,
job=2: draw confidence ellipse for prediction */
tf = 2. * (nobs-1) / (nobs * (nobs - 2));
ff = tf * FINV(prob,2,nobs-2);
if job = 2 then ff = ff * (nobs + 1.);
ff = sqrt(ff);
/* Eigenvectors are stored in columns */
CALL EIGEN(eval,evec,covm);
if eval[1] <= 0 then da = 0.;
else da = sqrt(eval[1]);
if eval[2] <= 0 then db = 0.;
else db = sqrt(eval[2]);
v0 = da * evec[1,1] * evec[1,1]
+ db * evec[1,2] * evec[1,2];
v1 = da * evec[1,1] * evec[2,1]
+ db * evec[1,2] * evec[2,2];
v2 = da * evec[2,1] * evec[2,1]
+ db * evec[2,2] * evec[2,2];
elp = j(nn+1,2,0.);
phi = 0.; dp = 1. / nn; s2 = 4. * sqrt(2.);
do j = 1 to nn+1;
if phi <= .25 then do;
v = phi - .125; xv = s2 * v;
yv = -sqrt(1. - xv * xv);
end; else
if phi <= .5 then do;
v = phi - .375; yv = s2 * v;
xv = sqrt(1. - yv * yv);
end; else
if phi <= .75 then do;
v = .625 - phi; xv = s2 * v;
yv = sqrt(1. - xv * xv);
end; else do;
v = .875 - phi; yv = s2 * v;
xv = -sqrt(1. - yv * yv);
end;
elp[j,1] = ff * (v0 * xv + v1 * yv) + xm[1];
elp[j,2] = ff * (v1 * xv + v2 * yv) + xm[2];
phi = phi + dp;
end;
/* Get axes of ellipsoid: a axis: 1 -> 2
b axis: 3 -> 4 */
if covm[1,1] > covm[2,2] then do;
t = da; da = db; db = t;
t = evec[1,1]; evec[1,1] = evec[1,2]; evec[1,2] = t;
t = evec[2,1]; evec[2,1] = evec[2,2]; evec[2,2] = t;
end;
rot = j(2,2,0.);
rot[1,1] = da * ff * evec[1,1];
rot[2,1] = da * ff * evec[1,2];
rot[1,2] = -db * ff * evec[1,2];
rot[2,2] = db * ff * evec[1,1];
corn = j(4,2,0.);
corn[1,1] = -1.; corn[1,2] = 0.;
corn[2,1] = 1.; corn[2,2] = 0.;
corn[3,1] = 0.; corn[3,2] = -1.;
corn[4,1] = 0.; corn[4,2] = 1.;
corn = corn * rot` + j(4,1,1) * xm;
finish elli2;
/* Initialization */
n = ncol(a); nobs = nrow(a);
nn = 120; /* number of ellipsoid points */
/* Compute classical scatter matrix */
muc = a[+,] / nobs;
covc = (a` * a - muc` * muc * nobs) / (nobs - 1.);
/* Compute robust scatter matrix */
CALL MVE(sc,xmve,dist,optn,a);
mur = xmve[1,]; covr = xmve[3:2+n,];
wgts = dist[3,];
/*--- do plots for all pairs of variables ---*/
ig = 0;
do jy = 2 to n;
do jx = 1 to jy-1;
ig = ig + 1;
call goptmve(ipsf,filn,n,jx,jy);
x = a[,jx]; y = a[,jy];
mu1 = (muc[jx] || muc[jy]);
cov1 = ( (covc[jx,jx] || covc[jx,jy]) //
(covc[jy,jx] || covc[jy,jy]) );
/* Compute points of classic Ellipsoid */
CALL ELLI1(job,prob,nobs,mu1,cov1,nn,corn1,elp1);
mu2 = (mur[jx] || mur[jy]);
cov2 = ( (covr[jx,jx] || covr[jx,jy]) //
(covr[jy,jx] || covr[jy,jy]) );
/* Compute points of robust Ellipsoid */
CALL ELLI1(job,prob,nobs,mu2,cov2,nn,corn2,elp2);
/* Get dimension of axes */
/* get [xmin,xmax] and [ymin,ymax] */
xmin1 = x[><]; xmax1 = x[<>];
ymin1 = y[><]; ymax1 = y[<>];
xmin2 = elp1[><,1]; xmax2 = elp1[<>,1];
ymin2 = elp1[><,2]; ymax2 = elp1[<>,2];
if xmin2 < xmin1 then xmin1 = xmin2;
if xmax2 > xmax1 then xmax1 = xmax2;
if ymin2 < ymin1 then ymin1 = ymin2;
if ymax2 > ymax1 then ymax1 = ymax2;
xmin2 = elp2[><,1]; xmax2 = elp2[<>,1];
ymin2 = elp2[><,2]; ymax2 = elp2[<>,2];
if xmin2 < xmin1 then xmin1 = xmin2;
if xmax2 > xmax1 then xmax1 = xmax2;
if ymin2 < ymin1 then ymin1 = ymin2;
if ymax2 > ymax1 then ymax1 = ymax2;
xd1 = xmax1 - xmin1; xc1 = .1 * xd1;
xmin2 = xmin1 - xc1; xmax2 = xmax1 + xc1;
xd2 = xmax2 - xmin2; xc2 = .1 * xd2;
xme = .5*(xmin2+xmax2);
yd1 = ymax1 - ymin1; yc1 = .1 * yd1;
ymin2 = ymin1 - yc1; ymax2 = ymax1 + yc1;
yd2 = ymax2 - ymin2; yc2 = .1 * yd2;
yme = .5*(ymin2+ymax2);
/* prepare the plot */
xbox = { 0 100 100 0};
ybox = { 0 0 100 100};
wind = (xmin2 || ymin2) // (xmax2 || ymax2);
/* create the plot */
call gstart;
call gopen("mve");
call gset("font","swiss"); /* set character font */
call gpoly(xbox,ybox); /* draw a box around plot */
call gwindow(wind); /* define window */
call gport({15 15 , 85 85}); /* define viewport */
call gxaxis((xmin2 || ymin2),xd2,11, , ,"4.1",1);
call gyaxis((xmin2 || ymin2),yd2,11, , ,"4.1",1);
/* plot xnam, ynam variable names of axis and the title */
call gset("height",1.5);
call gstrlen(len,vnam[jx]);
call gtext(xmax2-len,ymin2-yc2,vnam[jx]);
call gvtext(xmin2-xc2,ymax2,vnam[jy]);
call gset("height",2.0);
call gscenter(xme,ymax2,titl);
/* draw location point */
call gpoint(mu1[1],mu1[2],"circle","red",2);
/* draw points of ellipsoid */
call gpoint(elp1[,1],elp1[,2],"diamond","red",0.5);
/* plot the two ellipsoid axes */
/*
call gpoint(corn1[1,1],corn1[1,2],"diamond","red",1.2);
call gpoint(corn1[2,1],corn1[2,2],"diamond","red",1.2);
call gpoint(corn1[3,1],corn1[3,2],"diamond","red",1.2);
call gpoint(corn1[4,1],corn1[4,2],"diamond","red",1.2);
*/
call gdrawl(corn1[1,1:2],corn1[2,1:2], ,"red");
call gdrawl(corn1[3,1:2],corn1[4,1:2], ,"red");
/* draw location point */
call gpoint(mu2[1],mu2[2],"circle","green",2);
/* draw points of ellipsoid */
/* call gpoint(elp2[,1],elp2[,2],"dot","green",0.5); */
call gdraw(elp2[,1],elp2[,2], ,"green");
/* plot the two ellipsoid axes */
/*
call gpoint(corn2[1,1],corn2[1,2],"dot","green",1.2);
call gpoint(corn2[2,1],corn2[2,2],"dot","green",1.2);
call gpoint(corn2[3,1],corn2[3,2],"dot","green",1.2);
call gpoint(corn2[4,1],corn2[4,2],"dot","green",1.2);
*/
call gdrawl(corn2[1,1:2],corn2[2,1:2], ,"green");
call gdrawl(corn2[3,1:2],corn2[4,1:2], ,"green");
/* print the scatter points */
/* enumerate the zero weight points */
do i=1 to nobs;
xi = x[i]; yi = y[i];
if wgts[i] = 0 then do;
z = char(i,2,0);
call gpoint(xi,yi,"dot","red",1.5);
xi = xi + .01 * xd2;
call gscript(xi,yi,z,,,1.);
end;
else do;
call gpoint(xi,yi,"dot","green",1.5);
end; end;
call gshow;
call gclose;
call gstop;
end; end;
leave:
finish scatmve;
/*------------------------------------------------------------*/
start lmsdiap(job,optlms,optmve,prob,x,y,titl,ipsf,filn);
/* job = 1: plot LMS residuals vs. robust distances
= 2: plot LTS residuals vs. robust distances
= 15: do plot job=1 and job=5
= 25: do plot job=2 and job=5
ipsf: specifies GOPTIONS statement:
ipsf=0: draw interactively
ipsf=1: write file with smaller picture for document
ipsf=2: write file with larger picture for color printer
filn: prefix name of postscript file */
start gscenter(x,y,str);
call gstrlen(len,str);
call gscript(x-len/2,y,str);
finish gscenter;
start goptrd(job,ipsf,filn);
if ipsf then do;
if job=1 then post="lms";
if job=2 then post="lts";
if job=3 then post="ls";
fnam = (filn || post || ".ps");
fnam = rowcatc(fnam);
print fnam;
call execute("filename gsasfile'",fnam,"';");
if ipsf = 1 then do;
call execute('goptions reset=goptions dev=pslepsf
gsfmode=replace gaccess=gsasfile gsfname=GSASFILE
hsize=5.625 in vsize=4.65 in htext=3.0 pct border;');
end; else do;
call execute('goptions reset=goptions dev=pscolor
gsfmode=replace gaccess=gsasfile border erase;');
end;
end; else do;
call execute('goptions reset=goptions /*dev=xcolor*/;');
end;
finish goptrd;
start doplot(job,n,nobs,xdis,yres,wgts,titl);
if job=1 | job=2 then xstr= "Robust Distance RD";
if job=3 then xstr= "Classical Distance MD";
if job=1 then ystr= "LMS Std. Res.";
if job=2 then ystr= "LTS Std. Res.";
if job=3 then ystr= "LS Std. Res.";
chis = sqrt(cinv(.975,n)); /* print chis; */
xmin1 = xdis[><]; xmax1 = xdis[<>];
ymin1 = yres[><]; ymax1 = yres[<>];
if xmax1 < 1.1*chis then xmax1 = 1.1*chis;
if ymax1 < 3. then ymax1 = 3.;
if ymin1 > -3. then ymin1 = -3.;
xd1 = xmax1 - xmin1; xc1 = .1 * xd1;
xmin2 = xmin1 - xc1; xmax2 = xmax1 + xc1;
xd2 = xmax2 - xmin2; xc2 = .1 * xd2;
xme = .5*(xmin2+xmax2);
yd1 = ymax1 - ymin1; yc1 = .1 * yd1;
ymin2 = ymin1 - yc1; ymax2 = ymax1 + yc1;
yd2 = ymax2 - ymin2; yc2 = .1 * yd2;
yme = .5*(ymin2+ymax2);
/* prepare the plot */
xbox = { 0 100 100 0};
ybox = { 0 0 100 100};
wind = (xmin2 || ymin2) // (xmax2 || ymax2);
/* create the plot */
call gstart;
call gopen("diag");
call gset("font","swiss"); /* set character font */
call gpoly(xbox,ybox); /* draw a box around plot */
call gwindow(wind); /* define window */
call gport({15 15 , 85 85}); /* define viewport */
call gxaxis((xmin2 || ymin2),xd2,11, , ,"4.1",1);
call gyaxis((xmin2 || ymin2),yd2,11, , ,"4.1",1);
/* plot names of axis and the title */
call gset("height",1.5);
call gstrlen(xlen,xstr);
call gtext(xmax2-xlen,ymin2-yc2,xstr);
call gvtext(xmin2-xc2,ymax2,ystr);
call gset("height",2.0);
call gscenter(xme,ymax2,titl);
/* print the scatter points */
call gpoint(xdis,yres,"dot","magenta",1.5);
/* enumerate the zero weight points */
do i=1 to nobs;
if wgts[i] = 0 then do;
z = char(i,2,0);
xi = xdis[i]; yi = yres[i];
call gpoint(xi,yi,"dot","blue",1.5);
xi = xi + .01 * xd2;
call gscript(xi,yi,z,,,1.);
end; end;
xya = j(3,2,0.); xyb = j(3,2,0.);
xya[1,1]= xmin2; xya[1,2]= 2.5;
xyb[1,1]= xmax2; xyb[1,2]= 2.5;
xya[2,1]= xmin2; xya[2,2]=-2.5;
xyb[2,1]= xmax2; xyb[2,2]=-2.5;
xya[3,1]= chis; xya[3,2]=ymin2;
xyb[3,1]= chis; xyb[3,2]=ymax2;
call gdrawl(xya,xyb, ,"black");
call gshow;
call gclose;
call gstop;
finish doplot;
nobs = nrow(x); n = ncol(x);
if nrow(y) ^= nobs then do;
print "Data vectors x and y not compatible";
goto leave;
end;
/* Compute LMS or LTS: residuals are in (wgt[2,]); */
/* Compute LS: residuals are in (wgt[4,]); */
if job=1 | job=15 then do;
iob = 1;
optlms[3]= 2; /* ilsq: LS and LMS */
CALL LMS(sc,coef,wgt,optlms,y,x);
end;
if job=2 | job=25 then do;
iob = 2;
optlms[3]= 2; /* ilsq: LS and LTS */
CALL LTS(sc,coef,wgt,optlms,y,x);
end;
/* Compute robust distances: dist[2,] */
/* Mahalanobis Distances: dist[1,] */
CALL MVE(sc,xmve,dist,optmve,x);
/* Plot LMS_Res vs Rob_Dist */
xdis = dist[2,]; yres = wgt[2,]; wgts = dist[3,];
call goptrd(iob,ipsf,filn);
call doplot(iob,n,nobs,xdis,yres,wgts,titl);
if job=15 | job=25 then do;
/* Plot LS_Res vs Mah_Dist */
xdis = dist[1,]; yres = wgt[4,]; wgts = dist[3,];
call goptrd(3,ipsf,filn);
call doplot(3,n,nobs,xdis,yres,wgts,titl);
end;
leave:
finish lmsdiap;
/*------------------------------------------------------------*/
start scatmcd(job,optn,prob,a,vnam,titl,ipsf,filn);
/* job=1: draw confidence ellipse for mean,
job=2: draw confidence ellipse for prediction
ipsf: specifies GOPTIONS statement:
ipsf=0: draw interactively
ipsf=1: write file with smaller picture for document
ipsf=2: write file with larger picture for color printer
filn: prefix name of postscript file */
start gscenter(x,y,str);
call gstrlen(len,str);
call gscript(x-len/2,y,str);
finish gscenter;
start goptmcd(ipsf,filn,n,jx,jy);
if ipsf then do;
if n < 10 then do;
post = jy * 10 + jx;
post = char(post,2,0);
end; else do;
post = jy * 100 + jx;
post = char(post,4,0);
end;
fnam = (filn || post || ".ps");
fnam = rowcatc(fnam);
print fnam;
call execute("filename gsasfile'",fnam,"';");
if ipsf = 1 then do;
call execute('goptions reset=goptions dev=pslepsf
gsfmode=replace gaccess=gsasfile gsfname=GSASFILE
hsize=5.625 in vsize=4.65 in htext=3.0 pct border;');
end; else do;
call execute('goptions reset=goptions dev=pscolor
gsfmode=replace gaccess=gsasfile border erase;');
end;
end; else do;
call execute('goptions reset=goptions /*dev=xcolor*/;');
end;
finish goptmcd;
start elli1(job,prob,nobs,xm,covm,nn,corn,elp);
/* This is my version */
/* job=1: draw confidence ellipse for mean,
job=2: draw confidence ellipse for prediction */
tf = 2. * (nobs-1) / (nobs * (nobs - 2));
ff = tf * FINV(prob,2,nobs-2);
if job = 2 then ff = ff * (nobs + 1.);
ff = sqrt(ff);
/* Eigenvectors are stored in columns */
/* get orthogonal rotation matrix */
/* since length of eigenvectors is 1:
rot[1,1]= rot[2,2]= cos(alfa),
rot[2,1]= sin(alfa) = -rot[1,2] */
CALL EIGEN(eval,evec,covm);
if covm[1,1] > covm[2,2] then do;
t = eval[1]; eval[1] = eval[2]; eval[2] = t;
t = evec[1,1]; evec[1,1] = evec[1,2]; evec[1,2] = t;
t = evec[2,1]; evec[2,1] = evec[2,2]; evec[2,2] = t;
end;
if eval[1] <= 0 then da = 0.;
else da = sqrt(eval[1]);
if eval[2] <= 0 then db = 0.;
else db = sqrt(eval[2]);
rot = j(2,2,0.);
rot[1,1] = da * ff * evec[1,1];
rot[2,1] = da * ff * evec[1,2];
rot[1,2] = -db * ff * evec[1,2];
rot[2,2] = db * ff * evec[1,1];
/* Get axes of ellipsoid: a axis: 1 -> 2
b axis: 3 -> 4 */
corn = j(4,2,0.);
corn[1,1] = -1.; corn[1,2] = 0.;
corn[2,1] = 1.; corn[2,2] = 0.;
corn[3,1] = 0.; corn[3,2] = -1.;
corn[4,1] = 0.; corn[4,2] = 1.;
corn = corn * rot` + j(4,1,1) * xm;
/* Points of ellipsoid: elp[nn+1,2] */
pi = 3.141592653589; pi2 = 2. * pi;
elp = j(nn+1,2,0.);
phi = 0.; dp = pi2 / nn;
do j = 1 to nn+1;
x = cos(phi); y = sin(phi);
elp[j,] = (x || y) * rot` + xm;
phi = phi + dp;
end;
finish elli1;
start elli2(job,prob,nobs,xm,covm,nn,corn,elp);
/* This is the INSIGHT version */
/* job=1: draw confidence ellipse for mean,
job=2: draw confidence ellipse for prediction */
tf = 2. * (nobs-1) / (nobs * (nobs - 2));
ff = tf * FINV(prob,2,nobs-2);
if job = 2 then ff = ff * (nobs + 1.);
ff = sqrt(ff);
/* Eigenvectors are stored in columns */
CALL EIGEN(eval,evec,covm);
if eval[1] <= 0 then da = 0.;
else da = sqrt(eval[1]);
if eval[2] <= 0 then db = 0.;
else db = sqrt(eval[2]);
v0 = da * evec[1,1] * evec[1,1]
+ db * evec[1,2] * evec[1,2];
v1 = da * evec[1,1] * evec[2,1]
+ db * evec[1,2] * evec[2,2];
v2 = da * evec[2,1] * evec[2,1]
+ db * evec[2,2] * evec[2,2];
elp = j(nn+1,2,0.);
phi = 0.; dp = 1. / nn; s2 = 4. * sqrt(2.);
do j = 1 to nn+1;
if phi <= .25 then do;
v = phi - .125; xv = s2 * v;
yv = -sqrt(1. - xv * xv);
end; else
if phi <= .5 then do;
v = phi - .375; yv = s2 * v;
xv = sqrt(1. - yv * yv);
end; else
if phi <= .75 then do;
v = .625 - phi; xv = s2 * v;
yv = sqrt(1. - xv * xv);
end; else do;
v = .875 - phi; yv = s2 * v;
xv = -sqrt(1. - yv * yv);
end;
elp[j,1] = ff * (v0 * xv + v1 * yv) + xm[1];
elp[j,2] = ff * (v1 * xv + v2 * yv) + xm[2];
phi = phi + dp;
end;
/* Get axes of ellipsoid: a axis: 1 -> 2
b axis: 3 -> 4 */
if covm[1,1] > covm[2,2] then do;
t = da; da = db; db = t;
t = evec[1,1]; evec[1,1] = evec[1,2]; evec[1,2] = t;
t = evec[2,1]; evec[2,1] = evec[2,2]; evec[2,2] = t;
end;
rot = j(2,2,0.);
rot[1,1] = da * ff * evec[1,1];
rot[2,1] = da * ff * evec[1,2];
rot[1,2] = -db * ff * evec[1,2];
rot[2,2] = db * ff * evec[1,1];
corn = j(4,2,0.);
corn[1,1] = -1.; corn[1,2] = 0.;
corn[2,1] = 1.; corn[2,2] = 0.;
corn[3,1] = 0.; corn[3,2] = -1.;
corn[4,1] = 0.; corn[4,2] = 1.;
corn = corn * rot` + j(4,1,1) * xm;
finish elli2;
/* Initialization */
n = ncol(a); nobs = nrow(a);
nn = 120; /* number of ellipsoid points */
/* Compute classical scatter matrix */
muc = a[+,] / nobs;
covc = (a` * a - muc` * muc * nobs) / (nobs - 1.);
/* Compute robust scatter matrix */
CALL MCD(sc,xmcd,dist,optn,a);
mur = xmcd[1,]; covr = xmcd[3:2+n,];
wgts = dist[3,];
/*--- do plots for all pairs of variables ---*/
ig = 0;
do jy = 2 to n;
do jx = 1 to jy-1;
ig = ig + 1;
call goptmcd(ipsf,filn,n,jx,jy);
x = a[,jx]; y = a[,jy];
mu1 = (muc[jx] || muc[jy]);
cov1 = ( (covc[jx,jx] || covc[jx,jy]) //
(covc[jy,jx] || covc[jy,jy]) );
/* Compute points of classic Ellipsoid */
CALL ELLI1(job,prob,nobs,mu1,cov1,nn,corn1,elp1);
mu2 = (mur[jx] || mur[jy]);
cov2 = ( (covr[jx,jx] || covr[jx,jy]) //
(covr[jy,jx] || covr[jy,jy]) );
/* Compute points of robust Ellipsoid */
CALL ELLI1(job,prob,nobs,mu2,cov2,nn,corn2,elp2);
/* Get dimension of axes */
/* get [xmin,xmax] and [ymin,ymax] */
xmin1 = x[><]; xmax1 = x[<>];
ymin1 = y[><]; ymax1 = y[<>];
xmin2 = elp1[><,1]; xmax2 = elp1[<>,1];
ymin2 = elp1[><,2]; ymax2 = elp1[<>,2];
if xmin2 < xmin1 then xmin1 = xmin2;
if xmax2 > xmax1 then xmax1 = xmax2;
if ymin2 < ymin1 then ymin1 = ymin2;
if ymax2 > ymax1 then ymax1 = ymax2;
xmin2 = elp2[><,1]; xmax2 = elp2[<>,1];
ymin2 = elp2[><,2]; ymax2 = elp2[<>,2];
if xmin2 < xmin1 then xmin1 = xmin2;
if xmax2 > xmax1 then xmax1 = xmax2;
if ymin2 < ymin1 then ymin1 = ymin2;
if ymax2 > ymax1 then ymax1 = ymax2;
xd1 = xmax1 - xmin1; xc1 = .1 * xd1;
xmin2 = xmin1 - xc1; xmax2 = xmax1 + xc1;
xd2 = xmax2 - xmin2; xc2 = .1 * xd2;
xme = .5*(xmin2+xmax2);
yd1 = ymax1 - ymin1; yc1 = .1 * yd1;
ymin2 = ymin1 - yc1; ymax2 = ymax1 + yc1;
yd2 = ymax2 - ymin2; yc2 = .1 * yd2;
yme = .5*(ymin2+ymax2);
/* prepare the plot */
xbox = { 0 100 100 0};
ybox = { 0 0 100 100};
wind = (xmin2 || ymin2) // (xmax2 || ymax2);
/* create the plot */
call gstart;
call gopen("mcd");
call gset("font","swiss"); /* set character font */
call gpoly(xbox,ybox); /* draw a box around plot */
call gwindow(wind); /* define window */
call gport({15 15 , 85 85}); /* define viewport */
call gxaxis((xmin2 || ymin2),xd2,11, , ,"4.1",1);
call gyaxis((xmin2 || ymin2),yd2,11, , ,"4.1",1);
/* plot xnam, ynam variable names of axis and the title */
call gset("height",1.5);
call gstrlen(len,vnam[jx]);
call gtext(xmax2-len,ymin2-yc2,vnam[jx]);
call gvtext(xmin2-xc2,ymax2,vnam[jy]);
call gset("height",2.0);
call gscenter(xme,ymax2,titl);
/* draw location point */
call gpoint(mu1[1],mu1[2],"circle","red",2);
/* draw points of ellipsoid */
call gpoint(elp1[,1],elp1[,2],"diamond","red",0.5);
/* plot the two ellipsoid axes */
/*
call gpoint(corn1[1,1],corn1[1,2],"diamond","red",1.2);
call gpoint(corn1[2,1],corn1[2,2],"diamond","red",1.2);
call gpoint(corn1[3,1],corn1[3,2],"diamond","red",1.2);
call gpoint(corn1[4,1],corn1[4,2],"diamond","red",1.2);
*/
call gdrawl(corn1[1,1:2],corn1[2,1:2], ,"red");
call gdrawl(corn1[3,1:2],corn1[4,1:2], ,"red");
/* draw location point */
call gpoint(mu2[1],mu2[2],"circle","green",2);
/* draw points of ellipsoid */
/* call gpoint(elp2[,1],elp2[,2],"dot","green",0.5); */
call gdraw(elp2[,1],elp2[,2], ,"green");
/* plot the two ellipsoid axes */
/*
call gpoint(corn2[1,1],corn2[1,2],"dot","green",1.2);
call gpoint(corn2[2,1],corn2[2,2],"dot","green",1.2);
call gpoint(corn2[3,1],corn2[3,2],"dot","green",1.2);
call gpoint(corn2[4,1],corn2[4,2],"dot","green",1.2);
*/
call gdrawl(corn2[1,1:2],corn2[2,1:2], ,"green");
call gdrawl(corn2[3,1:2],corn2[4,1:2], ,"green");
/* print the scatter points */
/* enumerate the zero weight points */
do i=1 to nobs;
xi = x[i]; yi = y[i];
if wgts[i] = 0 then do;
z = char(i,2,0);
call gpoint(xi,yi,"dot","red",1.5);
xi = xi + .01 * xd2;
call gscript(xi,yi,z,,,1.);
end;
else do;
call gpoint(xi,yi,"dot","green",1.5);
end; end;
call gshow;
call gclose;
call gstop;
end; end;
leave:
finish scatmcd;
