/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 20.6                                              */
/*   TITLE: Figure 20.6 - Sequential T-optimum design                 */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "Figure 20.6 - Sequential T-optimum design";
%macro fig2006(sigma,seed=1,nAdd=41);

options nonotes;
   data Can;
      do x = -1 to 1 by 0.05;
         x10 = 1; x11 = exp( x); x12 = exp(-x);
         x20 = 1; x21 = x;       x22 = x*x;
         output;
         end;
   run;

   /*
   /  Initial design
   /------------------------------------------------------------------*/
   data Design;
      do x = -1 to 1;
         x10 = 1; x11 = exp( x); x12 = exp(-x);
         x20 = 1; x21 = x;       x22 = x*x;
         output;
         end;
   run;

   %let t10 =  4.5;
   %let t11 = -1.5;
   %let t12 = -2.0;

   data Eff;
      N = 1; x = -1; Eff = 0; output;
      N = 2; x =  0; Eff = 0; output;
   run;
   data _null_; set Design;
      call symput("NextX",trim(left(x)));
   run;
   %do i = 3 %to &nAdd;

      /*
      /  Compute Psi2 for current design.
      /---------------------------------------------------------------*/
      data Psi2; set Design;
         eta1 = &t10*x10 + &t11*x11 + &t12*x12;
      proc reg data=Psi2 noprint;
         model eta1 = x21 x22;
         output out=Psi2 r=r;
      data Psi2; set Psi2; Psi2 = r*r;
      proc summary data=Psi2;
         var psi2;
         output out=Del2 mean=Del2;
      data _null_; set Del2;
         call symput('Del2',trim(left(Del2)));
         call symput('Eff',trim(left(round(Del2/0.0010867,0.001))));
      run;
      data This;
         N = &i;
         x = &NextX;
         Eff = &Eff;
      run;
      %put &i: x=&NextX, Eff=&Eff;
	  data Eff; set Eff This;
      run;

      /*
      /  Simulate a new observation.
      /---------------------------------------------------------------*/
      data Design; set Design;
         Y = (&t10*x10 + &t11*x11 + &t12*x12) + &sigma*rannor(&seed);
      run;

      /*
      /  And add the point which maximizes Psi2 to the design.
      /---------------------------------------------------------------*/
      data Reg; set Design Can;
      proc reg data=Reg noprint; model y = x11 x12; output out=P1 p=Eta1;
      proc reg data=Reg noprint; model y = x21 x22; output out=P2 p=Eta2;
      data Psi2; merge P1 P2;
         where (y = .);
         Psi2 = (Eta1 - Eta2)**2;
      proc sort data=Psi2;
         by descending Psi2;
      data NextPoint; set Psi2(obs=1);
         call symput('NextX',trim(left(x)));
      data Design; set Design NextPoint(drop=Eta1 Eta2 Psi2);
      run;
      %end;

%mend;

title2 "Sigma = 0";
%fig2006(0,seed=1);
%let pagesize = %sysfunc(getoption(pagesize));
options ps=40;
proc plot data=Eff;
   plot Eff*N;
run;
options ps=&pagesize;
title2;

title2 "Sigma = 0.5";
%fig2006(0.5,seed=1);
%let pagesize = %sysfunc(getoption(pagesize));
options ps=40;
proc plot data=Eff;
   plot Eff*N;
run;
options ps=&pagesize;
title2;

title2 "Sigma = 1";
%fig2006(1,seed=1);
%let pagesize = %sysfunc(getoption(pagesize));
options ps=40;
proc plot data=Eff;
   plot Eff*N;
run;
options ps=&pagesize;
title2;

title2 "Sigma = 2";
%fig2006(2,seed=1);
%let pagesize = %sysfunc(getoption(pagesize));
options ps=40;
proc plot data=Eff;
   plot Eff*N;
run;
options ps=&pagesize;
title2;

title1;