/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 9.4                                               */
/*   TITLE: Optimum design measure for quadratic design on [-1,1]     */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*   NOTES: In order to modify this program to study designs for      */
/*          higher-order polynomials, add corresponding columns to    */
/*          the Zg matrix over the predictive grid, and perhaps also  */
/*          decrease the increment in the definition of the Grid      */
/*          data set.                                                 */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "Optimum design measure for quadratic design on [-1,1]";

/*
/  The approach below is discussed at length in Section 9.2 of the
/  text.
/---------------------------------------------------------------------*/
data Grid;
   do x = -1 to 1 by .1;
      output;
      end;
run;

proc iml;
   use Grid;
   read all var ("x") into xg;
   Ng = nrow(xg);
   Zg = j(Ng,1) || xg || xg##2;

   start WTrans(z);
      z2 = z##2;
      return(z2/sum(z2));
   finish;

   start DCritW(z) global(Zg);
      w = WTrans(z);
      return(log(det(Zg`*diag(w)*Zg)));
   finish;

   wInit = j(nrow(Zg),1) / nrow(Zg);
   call nlpqn(rc,        /* Return code                 */
              wOpt,      /* Returned optimal parameters */
              "DCritW",  /* Function to optimize        */
              wInit,     /* Initial value of parameters */
              1);        /* Specify a maximization      */
   wOpt = WTrans(wOpt)`;

   print "Raw optimized weights",,
         (xg || wOpt) [colname={"x" "w"}];

   print "Screening out tiny weights",,
         ((xg || wOpt)[loc(wOpt > 1e-4),]) [colname={"x" "w"}];
title1;