/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 6.1                                               */
/*   TITLE: Computations for Figures 6.1 & 6.2 - Confidence ellipses  */
/*    KEYS:                                                           */
/*    DATA: First six designs in Table 6.1                            */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "Computations for Figures 6.1 & 6.2 - Confidence ellipses";

/*
/  Create a single data set containing all designs.
/---------------------------------------------------------------------*/
data Table6_1;
   input Design nTrial @@;
   do i=1 to nTrial;
      input x @@;
      output;
      end;
   cards;
6.1 3 -1  0    1
6.2 6 -1 -1    0    0     1    1
6.3 8 -1 -1   -1   -1    -1   -1    1 1
6.4 5 -1 -0.5  0    0.5   1
6.5 7 -1 -1   -0.9 -0.85 -0.8 -0.75 1
6.6 2 -1  1
6.7 4 -1 -1    0    1
6.8 4 -1  0    0    1
;

proc iml;
   use Table6_1; read all var {Design x};
   free e;

/*
/  For each design ...
/---------------------------------------------------------------------*/
   do i = 1 to 6;

/*
/  ... create its variance matrix V ...
/---------------------------------------------------------------------*/
      l = loc(design(Design)[,i]);
      F = j(ncol(l),1) || x[l,];
      V = inv(F`*F);

/*
/  ... and compute the ellipse information from the eigenstructure
/  of M.  The major and minor axis lengths are the first two
/  eigenvalues, respectively; and the slope of the major axis is the
/  ratio of the elements of the first eigenvector.  Set the slope to
/  a missing value (.) if it's infinite.
/---------------------------------------------------------------------*/
      eval = eigval(V)`;
      evec = eigvec(V);
      if (evec[1,1]) then s = evec[2,1]/evec[1,1];
      else                s = .;
      e = e // ((6 + i/10) || eval || s);
      end;

/*
/  Print the results.
/---------------------------------------------------------------------*/
   print (1*e) [format=best6.4
                colname={"Design" "MajorAxis" "MinorAxis" "Slope"}];
quit;
title1;