/**********************************************************************/ /* Atkinson, Donev, and Tobias, "Optimum Experimental Designs" */ /* */ /* NAME: Program 20.1 */ /* TITLE: Example 20.1 & SAS Task 20.1 - Checking for second-order */ /* departures from a first-order model */ /* KEYS: */ /* DATA: */ /* NOTES: This program generates optimum continuous designs by */ /* optimizing only weights for given support, comparing */ /* results to optimum replications for large exact designs */ /* over the same support. */ /* */ /* Author: Randy Tobias */ /* History: */ /* Created...............................................27Jun2007 */ /**********************************************************************/ title1 "Example 20.1 - Checking for second-order departures from a first-order model"; %macro N0(alpha,N); /* Inverting equation 20.5: */ %if (%sysevalf(&alpha > 0)) %then %do; /* prior N0 corresponding */ %sysevalf(((1-&alpha)/&alpha)*&N) /* to given alpha */ %end; %else %do; %sysevalf(1000*&N) %end; %mend; data Grid; Source = "Grid"; do x1 = -1 to 1 by 0.05; do x2 = -1 to 1 by 0.05; x1 = round(x1,0.05); x2 = round(x2,0.05); output; end; end; run; /* / Use OPTEX to generate large exact designs for given alpha. /---------------------------------------------------------------------*/ %macro Exact(N,alpha,ds,wgt); proc optex data=Grid coding=none noprint; model x1 x2, x1*x1 x1*x2 x2*x2 / prior = 0,%N0((&alpha),&N); generate n=&N method=m_fedorov niter=10; output out=Design; proc freq data=Design noprint; table x2*x1 / out=&ds(rename=(Percent=&wgt)); run; %mend; %Exact(100,1/10 ,fDesign1,Wgt1); /* 100-run designs */ %Exact(100,2/ 3 ,fDesign2,Wgt2); %Exact(100,8/11 ,fDesign3,Wgt3); %Exact(100,0.7188,fDesign4,Wgt4); data fDesign100; merge fDesign1 fDesign2 fDesign3 fDesign4; by x2 x1; Wgt1 = Wgt1 / 100; label Wgt1 = "Alpha=0.1000"; Wgt2 = Wgt2 / 100; label Wgt2 = "Alpha=0.6667"; Wgt3 = Wgt3 / 100; label Wgt3 = "Alpha=0.7273"; Wgt4 = Wgt4 / 100; label Wgt4 = "Alpha=0.7188"; run; %Exact(1000,1/10 ,fDesign1,Wgt1); /* 1000-run designs */ %Exact(1000,2/ 3 ,fDesign2,Wgt2); %Exact(1000,8/11 ,fDesign3,Wgt3); %Exact(1000,0.7188,fDesign4,Wgt4); data fDesign1000; merge fDesign1 fDesign2 fDesign3 fDesign4; by x2 x1; Wgt1 = Wgt1 / 100; label Wgt1 = "Alpha=0.1000"; Wgt2 = Wgt2 / 100; label Wgt2 = "Alpha=0.6667"; Wgt3 = Wgt3 / 100; label Wgt3 = "Alpha=0.7273"; Wgt4 = Wgt4 / 100; label Wgt4 = "Alpha=0.7188"; run; %Exact(10000,1/10 ,fDesign1,Wgt1); /* 10000-run designs */ %Exact(10000,2/ 3 ,fDesign2,Wgt2); %Exact(10000,8/11 ,fDesign3,Wgt3); %Exact(10000,0.7188,fDesign4,Wgt4); data fDesign10000; merge fDesign1 fDesign2 fDesign3 fDesign4; by x2 x1; Wgt1 = Wgt1 / 100; label Wgt1 = "Alpha=0.1000"; Wgt2 = Wgt2 / 100; label Wgt2 = "Alpha=0.6667"; Wgt3 = Wgt3 / 100; label Wgt3 = "Alpha=0.7273"; Wgt4 = Wgt4 / 100; label Wgt4 = "Alpha=0.7188"; run; data Can; do x1 = -1 to 1 by 0.25; do x2 = -1 to 1 by 0.25; output; end; end; run; /* / Optimal weights for given alpha. /---------------------------------------------------------------------*/ proc iml; start MakeZ(x); Z = j(nrow(x),1) || x; do i = 1 to ncol(x); do j = i to ncol(x); Z = Z || x[,i]#x[,j]; end; end; return(Z); finish; start Crit(w,x) global(alpha); Z = MakeZ(x); K = diag((1:6)>3); Mi = Z`*diag(w)*Z; return(log(det( (1-alpha)*K /* Equation 20.4 and ff */ + alpha *Mi))); finish; start CritW(w) global(x); /* Crit for sqrt(w), x fixed */ return(Crit((w##2)/sum(w##2),x)); finish; /* / Quick and dirty, knowing that the support of the optimum design is / contained in the candidate set. /---------------------------------------------------------------------*/ start AlphaOpt; use Can; read all var {x1 x2} into xCan; x = xCan; /* Optimize weights for */ winit = sqrt(j(nrow(x),1) / nrow(x)); /* all support points. */ call nlpqn(rc,w,"CritW",winit,1); w = (w##2)`/sum(w##2); xOp = x[loc(w > 1e-5),]; /* Trim support points */ wOp = w[loc(w > 1e-5) ]; /* with negligible wgts. */ wOp = wOp /sum(wOp); x = xOp; /* And optimize again. */ winit = sqrt(wOp); call nlpqn(rc,w,"CritW",winit,1); w = (w##2)`/sum(w##2); xOp = x[loc(w > 1e-5),]; wOp = w[loc(w > 1e-5) ]; wOp = wOp /sum(wOp); wxOp = wOp || xOp; finish; alpha = 1/10; call AlphaOpt; create wDesign1 var {Wgt1 x1 x2}; append from wxOp; alpha = 2/ 3; call AlphaOpt; create wDesign2 var {Wgt2 x1 x2}; append from wxOp; alpha = 8/11; call AlphaOpt; create wDesign3 var {Wgt3 x1 x2}; append from wxOp; alpha = 0.7188; call AlphaOpt; create wDesign4 var {Wgt4 x1 x2}; append from wxOp; proc sort data=wDesign1; by x2 x1; proc sort data=wDesign2; by x2 x1; proc sort data=wDesign3; by x2 x1; proc sort data=wDesign4; by x2 x1; data wDesign; merge wDesign1 wDesign2 wDesign3 wDesign4; by x2 x1; label Wgt1 = "Alpha=0.1000"; label Wgt2 = "Alpha=0.6667"; label Wgt3 = "Alpha=0.7273"; label Wgt4 = "Alpha=0.7188"; run; title2 "Optimal frequencies for 100 runs"; proc print data=fDesign100 label noobs; var x2 x1 Wgt1 Wgt2 Wgt4 Wgt3; format Wgt1-Wgt4 6.4; run; title2; title2 "Optimal frequencies for 1000 runs"; proc print data=fDesign1000 label noobs; var x2 x1 Wgt1 Wgt2 Wgt4 Wgt3; format Wgt1-Wgt4 6.4; run; title2; title2 "Optimal frequencies for 10000 runs"; proc print data=fDesign10000 label noobs; var x2 x1 Wgt1 Wgt2 Wgt4 Wgt3; format Wgt1-Wgt4 6.4; run; title2; title2 "Optimal weights"; proc print data=wDesign label noobs; var x2 x1 Wgt1 Wgt2 Wgt4 Wgt3; format Wgt1-Wgt4 6.4; run; title2; title1;