/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 5.1                                               */
/*   TITLE: Table 5.4 - Determinants and variances for designs for a  */
/*          single quantitative factor                                */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "Table 5.4 - Determinants and variances for designs for a single quantitative factor";

/*
/  Define the three designs.
/---------------------------------------------------------------------*/
data Design5_1; input x @@; output; datalines;
-1 0 1
;
data Design5_2; input x @@; output; datalines;
-1 1 1
;
data Design5_3; input x @@; x = x / 3; output; datalines;
-3 -1 1 3
;

/*
/  Use IML to evaluate determinants and to maximize variance.
/---------------------------------------------------------------------*/
options nonotes;
proc iml;

/*
/  Define first and second order standardized prediction variances as
/  functions of a single value x; need to pull in the variance matrix
/  as a global.  The functions are written to handle vector input.
/---------------------------------------------------------------------*/
   start d1(x) global(V);
      z = j(nrow(x),1) || x;
      free d;
      do i = 1 to nrow(z);
         d = d // z[i,]*V*z[i,]`;
         end;
      return(d);
   finish;

   start d2(x) global(V);
      z = j(nrow(x),1) || x || x#x;
      free d;
      do i = 1 to nrow(z);
         d = d // z[i,]*V*z[i,]`;
         end;
      return(d);
   finish;

/*
/  A function to evaluate the maximum standardized prediction variance
/  over [-1,1], given the order of the model o and the design matrix
/  F.
/---------------------------------------------------------------------*/
   start dmax(o,F) global(V);
      V = inv(F`*F/nrow(F));

/*
/  Find the value of x on a grid over the region which maximizes d().
/---------------------------------------------------------------------*/
      xRegion = -1 + (1 - -1)*(0:100)`/100;
      if (o = 1) then dRegion = d1(xRegion);
      else            dRegion = d2(xRegion);
      xdmax = xRegion[dRegion[<:>]];

/*
/  Use a general nonlinear optimizer NLPQN to try to improve the
/  maximum value.  The penultimate argument specifies that the
/  optimizer should seek a maximum, and the last one gives bounds for
/  the region.
/---------------------------------------------------------------------*/
      if (o = 1) then call nlpqn(rc,xdmax,"d1",xdmax,1,{-1,1});
      else            call nlpqn(rc,xdmax,"d2",xdmax,1,{-1,1});

      if (o = 1) then return(d1(xdmax));
      else            return(d2(xdmax));
   finish;

/*
/  Evaluate determinants and max prediction variances.
/---------------------------------------------------------------------*/
   use Design5_1; read all var {x};
   F = j(nrow(x),1) || x;
   print "Design 5.1, First order model" (nrow(x)) (det(F`*F)) (dmax(1,F));

   use Design5_2; read all var {x};
   F = j(nrow(x),1) || x;
   print "Design 5.2, First order model" (nrow(x)) (det(F`*F)) (dmax(1,F));
                                                                      
   use Design5_1; read all var {x};
   F = j(nrow(x),1) || x || x#x;
   print "Design 5.1, Quadratic model  " (nrow(x)) (det(F`*F)) (dmax(2,F));

   use Design5_3; read all var {x};
   F = j(nrow(x),1) || x || x#x;
   print "Design 5.3, Quadratic model  " (nrow(x)) (det(F`*F)) (dmax(2,F));
quit;

title1;