Resources

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;