/**********************************************************************/
/*    Atkinson, Donev, and Tobias, "Optimum Experimental Designs"     */
/*                                                                    */
/*    NAME: Program 21.4                                              */
/*   TITLE: Example 21.3: First-order growth model                    */
/*    KEYS:                                                           */
/*    DATA:                                                           */
/*                                                                    */
/*  Author: Randy Tobias                                              */
/* History:                                                           */
/*    Created...............................................27Jun2007 */
/**********************************************************************/

title1 "Example 21.3: First-order growth model";

%let theta = 0.1;
%let nInc = 201;
data Can;
   t = 0;
      x = t*exp(-&theta*t);
      output;
   do it = 1 to &nInc;
      logt = log(1) + (it-1)*(log(30)-log(1))/(&nInc-1);
      t = exp(logt);
      x = t*exp(-(&theta)*t);
      output;
      end;
run;

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

   start DCritWT(wt) global(k);
      w = shape(wt[1:2],2,1);
      t = shape(wt[3:4],2,1);
      F1 = t#exp(-&theta*t);
      F = F1 || F1#F1;
      p = 2;
      r = 1;
      s = p - r;
      Dw = diag((w##2)/sum(w##2));
      m11 = F1`*Dw*F1 + 1e-8*i(ncol(F1));
      M   = F `*Dw*F  + 1e-8*i(ncol(F ));
      return((k/r)*log(det(M11)) + ((1-k)/s)*(log(det(M)) - log(det(M11))));
   finish;

   free kAll wAll tAll;
   do ik = 0 to 100;
      k = ik/100;
      wtinit =    {1 1      1 15};
      con    =    {. .      0  0}
               // {. .     30 30};
      call nlpqn(rc,wt,"DCritWT",wtinit,1,con);
      wOp = shape(wt[1:2],2,1);
      tOp = shape(wt[3:4],2,1);
      wOp = ((wOp##2)/sum(wOp##2));
      if (ik = 0) then do;
         wOp0 = wOp;
         tOp0 = tOp;
         end;
      else if (ik = 100) then do;
         wOp1 = wOp;
         tOp1 = tOp;
         end;
      kAll = kAll // k;
      wAll = wAll // wOp`;
      tAll = tAll // tOp`;
      end;
   print kAll wAll tAll;

   start D1(w,t);
      F1 = t#exp(-&theta*t);
      F = F1 || F1#F1;
      p = 2;
      r = 1;
      s = p - r;
      Dw = diag(w);
      m11 = F1`*Dw*F1;
      if (det(M11)<=0) then
         d1 = -999;
      else
         d1 = log(det(m11))/r;
      return(d1);
   finish;

   start Ds(w,t);
      F1 = t#exp(-&theta*t);
      F = F1 || F1#F1;
      p = 2;
      r = 1;
      s = p - r;
      Dw = diag(w);
      m11 = F1`*Dw*F1;
      M   = F `*Dw*F ;
      if (det(M) <= 0) then
         ds = -999;
      else if (det(M11) <= 0) then
         ds = 999;
      else
         ds = (log(det(M)) - log(det(M11)))/s;
      return(ds);
   finish;

   start EffD(w,t) global(wOp1,tOp1);
      d   = D1(w   ,t   );
      dOp = D1(wOp1,tOp1);
      return(exp(d)/exp(dOp));
   finish;

   start EffDs(w,t) global(wOp0,tOp0);
      d   = Ds(w   ,t   );
      dOp = Ds(wOp0,tOp0);
      return(exp(d)/exp(dOp));
   finish;

   free Eff;
   do i = 1 to nrow(kAll);
      w = wAll[i,]`;
      t = tAll[i,]`;
      Eff = Eff // (EffD(w,t) || EffDs(w,t));
      end;

   start Othert(t) global(x);
      y = &theta*t*exp(1-&theta*t);
      return((y - x)**2);
   finish;

   free t3;
   do i = 1 to nrow(kAll);
      x = &theta*tAll[i,1]*exp(1-&theta*tAll[i,1]);
      tinit = 15;
      con = {10 .}`;
      call nlpqn(rc,t,"Othert",tinit,0,con);
      t3 = t3 // t;
      end;

   data = kAll || wAll || tAll || Eff || t3;
   create fig21NL var {k w1 w2 t1 t2 EffD EffDs t3};
   append from data;

title2 "Fig. 21.3 - Compound optimum designs for parameter estimation and model checking";
%let pagesize = %sysfunc(getoption(pagesize));
options ps=40;
proc plot data=fig21Nl;
   plot t1*k;
   plot w1*k;
   plot EffD*k='-' EffDs*k='+' / overlay;
run;
options ps=&pagesize;
title2;

title1;