/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 21.2                                              */
/*   TITLE: Table 21.2 - Linear or quadratic regression               */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/


title1 "Table 21.2 - Linear or quadratic regression";

data Can;
   do x = -1 to 1 by 0.1;
      output;
      end;
run;

proc iml;
   use Can;
   read all var {x} into xCan;

   start MakeZ(x);
      F1 = j(nrow(x),1) || x;
      F2 = x#x;
      F = F1 || F2;
      return(F);
   finish;

   start llog(x); if (x > 0) then return(log(x)); else return(-9999); finish;

/*
/ Equation 21.18 in the text.
/---------------------------------------------------------------------*/
   start DCritW(w) global(F,k);
      p = ncol(F);
      r = 2;
      s = p - r;
      F1 = F[,1:r];
      Dw = diag((w##2)/sum(w##2));
      m11 = F1`*Dw*F1;
      M   = F `*Dw*F ;
	  z = 0;
	  if (  k) then z = z + (   k /r)* llog(det(M11));
	  if (1-k) then z = z + ((1-k)/s)*(llog(det(M)) - llog(det(M11)));
      return(exp(z));
   finish;

   start Optimize;
      x = XCan; F = MakeZ(x);
      winit = j(nrow(X),1);
      call nlpqn(rc,wOp,"DCritW",winit,1);
      wOp = ((wOp##2)/sum(wOp##2))`;
      xOp = x  [loc(wOp > 1e-3),];
      wOp = wOp[loc(wOp > 1e-3),];

      x = xOp; F = MakeZ(x);
      winit = sqrt(wOp);
      call nlpqn(rc,wOp,"DCritW",winit,1);
      wOp = ((wOp##2)/sum(wOp##2))`;
      xOp = x  [loc(wOp > 1e-3),];
      wOp = wOp[loc(wOp > 1e-3),];
   finish;

   k = 0  ; call Optimize; w1 = wOp; x1 = xOp;
   k = 0.5; call Optimize; w2 = wOp; x2 = xOp;
   k = 2/3; call Optimize; w3 = wOp; x3 = xOp;
   k = 1  ; call Optimize; w4 = wOp; x4 = xOp;

   k = 0  ;
      F = MakeZ(x1); D11 = DCritW(sqrt(w1));
      F = MakeZ(x2); D12 = DCritW(sqrt(w2));
      F = MakeZ(x3); D13 = DCritW(sqrt(w3));
      F = MakeZ(x4); D14 = DCritW(sqrt(w4));
   k = 0.5;
      F = MakeZ(x1); D21 = DCritW(sqrt(w1));
      F = MakeZ(x2); D22 = DCritW(sqrt(w2));
      F = MakeZ(x3); D23 = DCritW(sqrt(w3));
      F = MakeZ(x4); D24 = DCritW(sqrt(w4));
   k = 2/3;
      F = MakeZ(x1); D31 = DCritW(sqrt(w1));
      F = MakeZ(x2); D32 = DCritW(sqrt(w2));
      F = MakeZ(x3); D33 = DCritW(sqrt(w3));
      F = MakeZ(x4); D34 = DCritW(sqrt(w4));
   k = 1  ;
      F = MakeZ(x1); D41 = DCritW(sqrt(w1));
      F = MakeZ(x2); D42 = DCritW(sqrt(w2));
      F = MakeZ(x3); D43 = DCritW(sqrt(w3));
      F = MakeZ(x4); D44 = DCritW(sqrt(w4));

   k     = 0       // 0.5     // 2/3     // 1      ;
   w     = 2*w1[1] // 2*w2[1] // 2*w3[1] // 2*w4[1];
   wStar = (2 - k)/(4 - 3*k);
   EffD  = D41/D44 // D42/D44 // D43/D44 // D44/D44;
   EffDs = D11/D11 // D12/D11 // D13/D11 // D14/D11;
   print "D- and Ds-optimum designs for components",,
         k     [format=6.3]
         w     [format=6.3]
         wStar [format=6.3]
         EffD  [format=6.3]
         EffDs [format=6.3];

title1;