/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 6.5                                               */
/*   TITLE: Computations for Figure 6.9 - Variance dispersion graph   */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "Computations for Figure 6.9 - Variance dispersion graph";


/*
/  Create a grid of points on the unit circle.
/---------------------------------------------------------------------*/
data Circle;
   do i = 1 to 1000;
      t = (i-1)*(2*constant('pi'))/(1000-1);
      x1 = sin(t);
      x2 = cos(t);
      output;
      end;
run;

/*
/  Create the designs---a 2*2 and a 3*3 factorial, respectively.
/---------------------------------------------------------------------*/
proc factex;
   factors x1-x2;
   output out=Design1;
proc factex;
   factors x1-x2 / nlev=3;
   output out=Design2;
run;

proc iml;
   use Circle;  read all var ("x1":"x2") into sph;

/*
/  Create the variance matrix for the 2*2 design.
/---------------------------------------------------------------------*/
   use Design1; read all var ("x1":"x2") into des;
   n = nrow(des);
   F = j(n,1) || des;
   V = inv(F`*F);

/*
/  For selected radii ...
/---------------------------------------------------------------------*/
   free vdsp;
   do radius = 0 to 2 by 0.2;

/*
/  ... create the coefficients of the first-order model over the
/  circular grid at that radius.
/---------------------------------------------------------------------*/
      rsph = radius*sph;
      s = j(nrow(sph),1) || rsph;

/*
/  Compute then standardized variance for each grid point ...
/---------------------------------------------------------------------*/
      free d; do i = 1 to nrow(s); d = d // n*s[i,]*V*s[i,]`; end;

/*
/  ... and collect minimum, average, and maximum values.
/---------------------------------------------------------------------*/
      vdsp = vdsp // (radius || d[><] || d[:] || d[<>]);
      end;

/*
/  Print the results.  Note that since the design is rotatable, the
/  values for the minimum, mean, and maximum are all the same.
/---------------------------------------------------------------------*/
title2 "2*2 design and first-order model";
   print (1*vdsp) [colname={"Radius" "VMin" "VAvg" "VMax"} format=8.2];
title2;



proc iml;
   use Circle;  read all var ("x1":"x2") into sph;

/*
/  Create the variance matrix for the 3*3 design.
/---------------------------------------------------------------------*/
   use Design2; read all var ("x1":"x2") into des;
   n = nrow(des);
   F = j(n,1) || des;
   do i = 1 to ncol(des);
      do j = i to ncol(des);
         F = F || des[,i]#des[,j];
         end;
      end;
   V = inv(F`*F);

/*
/  For selected radii ...
/---------------------------------------------------------------------*/
   free vdsp;
   do radius = 0 to 2 by 0.2;

/*
/  ... create the coefficients of the second-order model over the
/  circular grid at that radius.
/---------------------------------------------------------------------*/
      rsph = radius*sph;
      s = j(nrow(rsph),1) || rsph;
      do i = 1 to ncol(rsph);
         do j = i to ncol(rsph);
            s = s || rsph[,i]#rsph[,j];
            end;
         end;

/*
/  Compute then standardized variance for each grid point ...
/---------------------------------------------------------------------*/
      free d; do i = 1 to nrow(s); d = d // n*s[i,]*V*s[i,]`; end;

/*
/  ... and collect minimum, average, and maximum values.
/---------------------------------------------------------------------*/
      vdsp = vdsp // (radius || d[><] || d[:] || d[<>]);
      end;

/*
/  Print the results.
/---------------------------------------------------------------------*/
title2 "3*3 design and second-order model";
   print (1*vdsp) [colname={"Radius" "VMin" "VAvg" "VMax"} format=8.2];
title2;

title1;