/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 11.2                                              */
/*   TITLE: Example 11.1 - Cubic regression through the origin        */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "Figure 11.6 & 11.7 - Cubic regression through the origin";


/*
/  Create candidates [0,1] by .05.
/---------------------------------------------------------------------*/
data Candidates; do x = 0 to 1 by .05; output; end;

/*
/  Create initial 3-point design.
/---------------------------------------------------------------------*/
data Design3;
   x = 0.1; output;
   x = 0.5; output;
   x = 0.9; output;
run;

/*
/  In order to add points to a design, name the previous design with
/  the AUGMENT= option and give the new design size.
/---------------------------------------------------------------------*/
title2 "Adding the first point";
proc optex data=Candidates coding=none;
   model x x*x x*x*x / noint;
   generate n=4 augment=Design3 method=sequential;
   examine points;
run;
title2;

/*
/  To explore how sequential augmentation behaves as more points are
/  added, the following macro adds one point at a time and saves the
/  point that was added and the efficiencies of the augmented design.
/---------------------------------------------------------------------*/
%macro AddPoints;

options nonotes;

data Last; set Design3; run;

data Added; if (0); run;

%do N = 3 %to 50;
   data Candidates; set Candidates; Source = 'Can ';
   data Last;       set Last;       Source = 'Last';

   ods listing close;
   proc optex data=Candidates coding=none;
      model x x*x x*x*x / noint;
      generate initdesign=Last method=sequential;
      ods output Efficiencies=Eff;
   run;

   proc optex data=Candidates coding=none;
      model x x*x x*x*x / noint;
      generate n=%eval(&N+1) augment=Last method=sequential;
      output out=Next;
      id Source;
   quit;
   ods listing;

   data Add; merge Next(where=(Source ^= 'Last')) Eff;
      N = &N;
      drop Source DesignNumber;
   data Added; set Added Add;
   data Last; set Next;
   run;
   %end;

options notes;
%mend;

%AddPoints;

/*
/  Compute the D-criterion for the optimal continuous design.
/---------------------------------------------------------------------*/
proc iml;

/*
/  The algebraic optimum
/---------------------------------------------------------------------*/
   xO =    (5-sqrt(5))/10
        // (5+sqrt(5))/10
        // 1;
   ZO =    xO
        || xO#xO
        || xO#xO#xO;
   MOpt = ZO`*ZO/3;
   DOpt = det(MOpt)**(1/3);

/*
/  The numerical optimum on 100 equally-spaced values of x in [0,1]
/---------------------------------------------------------------------*/
   start DCritW(w) global(Z);
      Dw = diag((w##2)/sum(w##2));
      return(log(det(Z`*Dw*Z)**(1/ncol(Z))));
   finish;

   nintvl = 100;
   winit = j(nintvl,1) / nintvl;
   x = 0 + ((1:nintvl)` - 1)*(1 - 0)/(nintvl - 1);
   Z =    x
       || x#x
       || x#x#x;
   call nlpqn(rc,w,"DCritW",winit,1);

/*
/  Noting that the optimum weights seem to point to equal weighting
/  over a three-point support, assume three points and look for the
/  optimum three points.
/---------------------------------------------------------------------*/
   start DCritX(xx);
      x = shape(xx,3,1);
      Z =    x
          || x#x
          || x#x#x;
      d = det(Z`*Z/3)**(1/3);
      if (d < exp(-10)) then return(-10);
      else                   return(log(d));
   finish;

   title2 "Finding the exact D-optimal design";
   xinit = (1:3)/3;
   con   = {0 0 0,  /* All three points constrained to [0,1] */
            1 1 1};
   call nlpqn(rc,x,"DCritX",xinit,{1 1},con);

   print (x`) (exp(DCritX(x))) DOpt;
   title2;

   call symput('DOpt',trim(left(char(exp(DCritX(x)),20))));

/*
/  OPTEX does not make the intermediate pieces of the optimization
/  available, but for the purposes of this table they can be
/  reconstructed.  For example, we can use (?) to get d(Xi) from the
/  successive optimal log-determinants.
/---------------------------------------------------------------------*/
data dXi; set Added;
   dXi  = (N-1)*(exp(DCriterion - lag(DCriterion))-1);
   if (dXi ^= .);
   keep dXi;

/*
/  OPTEX's G-criterion is the square-root (expressed as a percentage)
/  of (?), and its log-determinant D-criterion needs to be scaled and
/  transformed appropriately to compare with the optimal D-criterion.
/---------------------------------------------------------------------*/
data Table11_1b; merge Added dXi;
   GEff = (GCriterion/100)**2;
   DEff = exp(DCriterion/3 - log(N))/&DOpt;
Proc print data=Table11_1b noobs;
   var N x dXi GEff DEff;
   format GEff DEff 5.3;
   format dXi 5.2;
   format x   best4.;
run;