/**********************************************************************/ /* Atkinson, Donev, and Tobias, "Optimum Experimental Designs" */ /* */ /* NAME: Program 19.5 */ /* TITLE: SAS Task 19.2 - Augmenting the design of Fig. 19.3 */ /* KEYS: */ /* DATA: */ /* */ /* Author: Randy Tobias */ /* History: */ /* Created...............................................27Jun2007 */ /**********************************************************************/ title1 "Table 19.1 - Augmenting the design of Fig. 19.3"; data Initial; input x1 x2; cards; -1 -1 -0.5 -1 0 -1 -1 0 0 0 0 -0.5 ; data Candidates; do x2 = -1 to 1 by 0.1; do x1 = -1 to 1 by 0.1; output; end; end; run; /* / Optimization with respect to the full 20x20 region takes more than / half an hour: it speeds things up considerably to consider only / points around -1,0,1. /---------------------------------------------------------------------*/ data Candidates; set Candidates; x1 = round(x1,0.1); x2 = round(x2,0.1); where (abs(x1 - round(x1)) <= 0.2+1e-6) & (abs(x2 - round(x2)) <= 0.2+1e-6); run; 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; /* / Read in intial and candidate points and create their respective / design matrices. /---------------------------------------------------------------------*/ use Initial; read all var {x1 x2} into x; n0 = nrow(x); FInit = MakeZ(x); use Candidates; read all var {x1 x2} into x; nCan = nrow(x); FCan = MakeZ(x); start DCritW(w) global(FInit,F,alpha,n0,nCan); Dw = diag((w##2)/sum(w##2)); return(log(det( (1-alpha)*FInit` *FInit/n0 + alpha *F `*Dw*F ))); finish; /* / A quick-and-dirty approach: Optimize twice, sweeping out tiny / weights. /---------------------------------------------------------------------*/ start AlphaOpt; F = FCan; winit = j(nCan,1) / nCan; call nlpqn(rc,w,"DCritW",winit,1); w = ((w##2)/sum(w##2))`; F = F[loc(w > 1e-5),]; winit = w[loc(w > 1e-5),]; call nlpqn(rc,w,"DCritW",winit,1); w = ((w##2)/sum(w##2))`; FOp = F[loc(w > 1e-5),]; wOp = w[loc(w > 1e-5),]; wOp = wOp/sum(wOp); finish; free FAll wAll Which; /* / Compute optimum for different alpha and collect results. /---------------------------------------------------------------------*/ alpha = 1/3; call AlphaOpt; Which = Which // 1*j(nrow(FOp),1) ; Fall = Fall // FOp[,2:3]; Wall = Wall // Wop ; alpha = 3/5; call AlphaOpt; Which = Which // 2*j(nrow(FOp),1) ; Fall = Fall // FOp[,2:3]; Wall = Wall // Wop ; alpha = 13/19; call AlphaOpt; Which = Which // 3*j(nrow(FOp),1) ; Fall = Fall // FOp[,2:3]; Wall = Wall // Wop ; alpha = 19/25; call AlphaOpt; Which = Which // 4*j(nrow(FOp),1) ; Fall = Fall // FOp[,2:3]; Wall = Wall // Wop ; alpha = 1; call AlphaOpt; Which = Which // 5*j(nrow(FOp),1) ; Fall = Fall // FOp[,2:3]; Wall = Wall // Wop ; data = Which || FAll || wAll; create FAll var {Which x1 x2 w}; append from data; quit; /* / Print collected results. /---------------------------------------------------------------------*/ data FAll; set FAll; format w 5.3; proc sort data=FAll; by Which x1 x2; data F1; set FAll; where (Which=1); drop Which; rename w=W1; R1 = round( 3*w); data F2; set FAll; where (Which=2); drop Which; rename w=W2; R2 = round( 9*w); data F3; set FAll; where (Which=3); drop Which; rename w=W3; R3 = round(13*w); data F4; set FAll; where (Which=4); drop Which; rename w=W4; R4 = round(19*w); data F5; set FAll; where (Which=5); drop Which; rename w=W5; R5 = round(13*w); data Table1901; merge F1 F2 F3 F4 F5; by x1 x2; proc sort data=Table1901; by x2 x1; data Table1901; set Table1901; Constant = 1; run; title2 "Continuous optimum designs and integer approximations"; proc print data=Table1901 noobs; var x1 x2 w1 r1 w2 r2 w3 r3 w4 r4 w5 r5; format w1-w5 5.3; run; title2; /* / For comparison, compute the true exact optimum designs / corresponding to each of the continuous optimum designs. /---------------------------------------------------------------------*/ data Candidates; set Candidates; Source = "Candidates"; data Prior; set Initial; Source = "Initial "; do i = 1 to 10; output; end; run; title2 "Optimum design for nAugment = 3"; data Prior; set Initial; Source = "Initial "; proc optex data=Candidates coding=orthcan; model x1 x2 x1*x1 x1*x2 x2*x2; generate n=9 augment=Prior method=m_fedorov niter=10000 keep=10; output out=Design; id Source; proc freq data=Design noprint; table Source*x1*x2 / list out=G3; run; title2; title2 "Optimum design for nAugment = 9"; proc optex data=Candidates coding=orthcan; model x1 x2 x1*x1 x1*x2 x2*x2; generate n=15 augment=Prior method=m_fedorov niter=10000 keep=10; output out=Design; id Source; proc freq data=Design noprint; table Source*x1*x2 / list out=G9; run; title2; title2 "Optimum design for nAugment = 13"; proc optex data=Candidates coding=orthcan; model x1 x2 x1*x1 x1*x2 x2*x2; generate n=19 augment=Prior method=m_fedorov niter=10000 keep=10; output out=Design; id Source; proc freq data=Design noprint; table Source*x1*x2 / list out=G13; run; title2; title2 "Optimum design for nAugment = 19"; proc optex data=Candidates coding=orthcan; model x1 x2 x1*x1 x1*x2 x2*x2; generate n=25 augment=Prior method=m_fedorov niter=10000 keep=10; output out=Design; id Source; proc freq data=Design noprint; table Source*x1*x2 / list out=G19; run; title2; title2 "Optimum design for nAugment = Infinity"; proc optex data=Candidates coding=orthcan; model x1 x2 x1*x1 x1*x2 x2*x2; generate n=13 method=m_fedorov niter=10000 keep=10; output out=Design; id Source; proc freq data=Design noprint; table Source*x1*x2 / list out=GOpt; run; title2; title2 "Optimum exact designs"; data G3 ; set G3 ; where (Source="Candidates"); drop Source PERCENT; rename COUNT=W3 ; data G9 ; set G9 ; where (Source="Candidates"); drop Source PERCENT; rename COUNT=W9 ; data G13 ; set G13 ; where (Source="Candidates"); drop Source PERCENT; rename COUNT=W13 ; data G19 ; set G19 ; where (Source="Candidates"); drop Source PERCENT; rename COUNT=W19 ; data GOpt; set GOpt; where (Source="Candidates"); drop Source PERCENT; rename COUNT=WOpt; data OP; merge G3 G9 G13 G19 GOpt; by x1 x2; proc sort data=Op; by x2 x1; proc print data=Op noobs; run; title2; title1;