/**********************************************************************/ /* Atkinson, Donev, and Tobias, "Optimum Experimental Designs" */ /* */ /* NAME: Program 12.3 */ /* TITLE: Example 12.1 - Exact optimal quadratic designs */ /* KEYS: */ /* DATA: */ /* NOTES: Same as Program 13.1 - SAS Task 13.1 */ /* */ /* Author: Randy Tobias */ /* History: */ /* Created...............................................27Jun2007 */ /**********************************************************************/ title1 "Example 12.1 continued - Exact optimal quadratic designs"; title2 "Using OPTEX to compute exact optimal designs"; proc factex; factors x1 x2 / nlev=3; output out=Grid; proc optex data=Grid coding=orthcan; model x1|x2 x1*x1 x2*x2; generate n=6 method=m_fedorov niter=1000 keep=10; examine points; output out=Grid6; title3 "N = 6"; run; generate n=7 method=m_fedorov niter=1000 keep=10; examine points; output out=Grid7; title3 "N = 7"; run; generate n=8 method=m_fedorov niter=1000 keep=10; examine points; output out=Grid8; title3 "N = 8"; run; generate n=9 method=m_fedorov niter=1000 keep=10; examine points; output out=Grid9; title3 "N = 9"; run; generate n=13 method=m_fedorov niter=1000 keep=10; examine points; output out=Grid13; title3 "N = 13"; run; title2; /* / Use constraied nonlinear optimization to find optimal designs / directly. Begin by defining a function to compute the log / determininant for a quadratic model. /---------------------------------------------------------------------*/ options nonotes; proc iml; start llog(x); if (x > 0) then return(log(x)); else return(-9999); finish; start Crit2(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; Mi = Z`*Z/nrow(x); return(llog(det(Mi))/ncol(Mi)); finish; start CritX(_x) global(m,N); return(Crit2(shape(_x,N,m))); finish; start tlc(n); return(trim(left(char(n,12)))); finish; /* / A function to find the best exact design for a given number of / factors and points. The approach uses 10 random restarts; this / might have to be modified for more than 2 factors. /---------------------------------------------------------------------*/ if (inputn(symget('sysver'),'3.1') >= 9) then do; start sort(x,idx); sortx = j(nrow(x), ncol(x), .); do i = nrow(idx)*ncol(idx) to 1 by -1; sortx[rank(x[,idx[i]]),] = x[]; end; x = sortx; finish; end; start OpDesign(_m,_N) global(N,m); N = _N; m = _m; BestD = -9999; free BestX; do i = 1 to 100; xinit = 2*ranuni(1:(N*m))-1; con = shape(-1,1,N*m) // shape( 1,1,N*m); call nlpqn(rc,xOp,"CritX",xinit,1,con); if (CritX(xOp) > BestD) then do; BestD = CritX(xOp); BestX = xOp; end; end; x = shape(BestX,N,m); call sort(x,1:ncol(x)); return(x); finish; /* / Read in grid designs and compute corresponding designs over the / continuous region. /---------------------------------------------------------------------*/ use Grid6; read all var {x1 x2} into Grid6; close; use Grid7; read all var {x1 x2} into Grid7; close; use Grid8; read all var {x1 x2} into Grid8; close; use Grid9; read all var {x1 x2} into Grid9; close; use Grid13; read all var {x1 x2} into Grid13; close; Cont6 = OpDesign(2,6); Cont7 = OpDesign(2,7); Cont8 = OpDesign(2,8); Cont9 = OpDesign(2,9); Cont13 = OpDesign(2,13); /* / Compute optimal design weights: see Program 12.2. /---------------------------------------------------------------------*/ start Crit(w,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; Mi = Z`*diag(w)*Z; return(log(det(Mi))/ncol(Mi)); finish; start CritW(w) global(x); return(Crit((w##2)/sum(w##2),x)); finish; start Consolidate(x,w,ex,ew); x = x[loc(w > ew),]; w = w[loc(w > ew)]; Group = shape(0,1,ncol(x)); do i = 1 to ncol(x); if (^Group[i]) then do; Group[i] = max(Group) + 1; do j = i+1 to ncol(x); if (max(abs(x[i,] - x[j,])) < ex) then Group[j] = Group[i]; end; end; end; xNew = shape(.,ncol(x),max(Group)); wNew = shape(.,1 ,max(Group)); do i = 1 to max(Group); xNew[i,] = (x[loc(Group = i),])[,:]; wNew[i ] = (w[loc(Group = i) ])[ +]; end; wx = w || x; call sort(wx,1:ncol(x)+1); w = wx[,1]; x = wx[,2:ncol(x)+1]; finish; start tlc(n); return(trim(left(char(n,12)))); finish; use Grid; read all var {"x1" "x2"} into xGrid; x = xGrid; winit = j(nrow(x),1) / nrow(x); winit = sqrt(winit); call nlpqn(rc,wOp,"CritW",winit,1); wOp = (wOp##2)`/sum(wOp##2); xOp = x; call Consolidate(xOp,wOp,0,1e-4); x = xOp; winit = sqrt(wOp); call nlpqn(rc,wOp,"CritW",winit,1); wOp = (wOp##2)`/sum(wOp##2); xOp = x; start DEff(x) global(wOp,xOp); w = j(nrow(x),1)/nrow(x); return(exp(Crit(w,x))/exp(Crit(wOp,xOp))); finish; N = ( 6 || 7 || 8 || 9 || 13 )`; Grid = (DEff(Grid6)||DEff(Grid7)||DEff(Grid8)||DEff(Grid9)||DEff(Grid13))`; Cont = (DEff(Cont6)||DEff(Cont7)||DEff(Cont8)||DEff(Cont9)||DEff(Cont13))`; Ratio = Grid/Cont; /* / The results confirm that Box and Draper's designs are D-optimal. /---------------------------------------------------------------------*/ print "Efficiencies of exact designs over 3x3 grid and [-1,1]^2",, N [format=best2.] Grid [format=7.4] Cont [format=7.4] Ratio [format=7.4]; /* / Here's a little trick to make OPTEX compute its efficiencies / with respect to the true optimum design weight. First, note that / multiplying the optimal design weights by 55 results in near / integer values. Thus, the design with N=sum(round(55*wOp)) runs / and replications given by round(55*wOp) will be nearly optimal. / In fact, the following calculation shows it to be fully efficient. /---------------------------------------------------------------------*/ iOp = round(55*wOp/wOp[1]); N = sum(iOp); free xiOp; do i = 1 to nrow(iOp); do j = 1 to iOp[i]; xiOp = xiOp // xOp[i,]; end; end; print "Creating an optimal replicated design",, N iOp (iOp/N) wOp (exp(Crit(iOp/N,xOp))/exp(Crit(wOp,xOp))) [format=percent9.2] (DEff(xiOp)) [format=percent9.2]; /* / Recall that with CODING=ORTHCAN, OPTEX's efficiencies are / effectively computed with respect to a design that's uniform over / the candidates. Thus, if you use this optimally replicated design / as the candidate set, the efficiencies that OPTEX computes will be / very close to those above. /---------------------------------------------------------------------*/ create iOp var {"x1" "x2"}; append from xiOp; title2 "Using OPTEX to compute true efficiencies"; proc optex data=iOp coding=orthcan; model x1|x2 x1*x1 x2*x2; generate initdesign=Grid6 method=sequential; ods select Efficiencies; title3 "N = 6"; run; generate initdesign=Grid7 method=sequential; ods select Efficiencies; title3 "N = 7"; run; generate initdesign=Grid8 method=sequential; ods select Efficiencies; title3 "N = 8"; run; generate initdesign=Grid9 method=sequential; ods select Efficiencies; title3 "N = 9"; run; generate initdesign=Grid13 method=sequential; ods select Efficiencies; title3 "N = 13"; run; title3; title2; title1;