/**********************************************************************/ /* Atkinson, Donev, and Tobias, "Optimum Experimental Designs" */ /* */ /* NAME: Program 18.1 */ /* TITLE: Example 18.1: Truncated quadratic model */ /* KEYS: */ /* DATA: */ /* */ /* Author: Randy Tobias */ /* History: */ /* Created...............................................27Jun2007 */ /**********************************************************************/ title1 "Example 18.1: Truncated quadratic model"; proc iml; /* / A little function to turn a number into a string for printing. /------------------------------------------------------------------*/ start tlc(n); return(trim(left(char(n,12)))); finish; /* / This function simplifies a design, often revealing an equivalent / design with fewer support and thus one whose structure is easier / to see. /------------------------------------------------------------------*/ start Consolidate(x,w,ex,ew); x = x[loc(w > ew)]`; /* Drop points with small wgts */ w = w[loc(w > ew)]`; Group = shape(0,1,ncol(x)); /* Replace nearby points with */ do i = 1 to ncol(x); /* their avg, with wgts equal */ if (^Group[i]) then do; /* to sum(wgt) over averaged */ Group[i] = max(Group) + 1; /* points. */ do j = i+1 to ncol(x); if (abs(x[i] - x[j]) < ex) then Group[j] = Group[i]; end; end; end; xNew = shape(.,1,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; r = rank(xNew); /* Make sure the points are */ x = xNew; /* sorted by x. */ w = wNew; x[r] = xNew; w[r] = wNew; finish; /* / Define 5-point prior as given in the text. /---------------------------------------------------------------------*/ tPop = {0.3 0.6 1 1.5 2}; pPop = j(1,ncol(tPop))/ncol(tPop); start llog(x); if (x > 0) then return(log(x)); else return(-9999); finish; /* / The approach is much the same as for optimum weighted designs in / previous chapters, except that the D-criterion is now an average. / The particular nature of the average is defined in Table 18.1 of / the text. /---------------------------------------------------------------------*/ start DCrit(w,x) global(tPop,pPop,criterion); if (criterion = 1) then EldMi = 0; else if (criterion = 2) then E_dMi = 0; else if (criterion = 3) then E__Mi = shape(0,1,1); else if (criterion = 4) then E_dM_ = 0; else E__M_ = shape(0,1,1); do i = 1 to ncol(tPop); fi = tPop[i]*(x#(1 - tPop[i]*x))#(0 < x)#(x < 1/tPop[i]); Mi = fi*diag(w)*fi`; if (criterion = 1) then EldMi = EldMi + pPop[i]*llog(det(ginv(Mi))); else if (criterion = 2) then E_dMi = E_dMi + pPop[i]* det(ginv(Mi)); else if (criterion = 3) then E__Mi = E__Mi + pPop[i]* ginv(Mi) ; else if (criterion = 4) then E_dM_ = E_dM_ + pPop[i]* det( Mi ); else E__M_ = E__M_ + pPop[i]* Mi ; end; if (criterion = 1) then return(- EldMi ); else if (criterion = 2) then return(-llog( E_dMi )); else if (criterion = 3) then return(-llog(det(E__Mi))); else if (criterion = 4) then return( llog( E_dM_ )); else return( llog(det(E__M_))); finish; start DCritW(w) global(x); return(DCrit((w##2)/sum(w##2),x)); finish; start DCritX(x) global(w); return(DCrit(w,x)); finish; start DCritWX(wx); w = shape(wx,2,ncol(wx)/2)[1,]; w = (w##2)/sum(w##2); x = shape(wx,2,ncol(wx)/2)[2,]; return(DCrit(w,x)); finish; criterion = 2; /* Criterion II: log E | inv(M) | */ /* / Compute optimum design weights over 10-point grid. /---------------------------------------------------------------------*/ nint = 10; x = (1:nint)/nint; winit = j(ncol(x),1) / ncol(x); call nlpqn(rc,w10,"DCritW",winit,1); w10 = ((w10##2)/sum(w10##2)); x10 = x; call Consolidate(x10,w10,0,1e-4); /* / Compute optimum design weights over 20-point grid. /---------------------------------------------------------------------*/ nint = 20; x = (1:nint)/nint; winit = j(ncol(x),1) / ncol(x); call nlpqn(rc,w20,"DCritW",winit,1); w20 = ((w20##2)/sum(w20##2))`; x20 = x; call Consolidate(x20,w20,0,1e-4); /* / Compute optimum design support and weights over [0,1], starting / from the 50-point grid. /---------------------------------------------------------------------*/ nint = 50; x = (1:nint)/nint; winit = j(ncol(x),1) / ncol(x); call nlpqn(rc,wGrid,"DCritW",winit,1); wGrid = ((wGrid##2)/sum(wGrid##2)); xGrid = x [loc(wGrid > 0)]; call Consolidate(xGrid,wGrid,0,1e-4); xOp = xGrid; wOp = wGrid; i = 0; Improvement = 1; do while ((i < 100) & (Improvement > 1e-4)); i = i + 1; last = (xOp` || wOp`); wxinit = sqrt(wOp) || xOp; con = (shape(.,1,ncol(xOp)) || shape(0,1,ncol(xOp))) // (shape(.,1,ncol(xOp)) || shape(1,1,ncol(xOp))); call nlpqn(rc,wxOp,"DCritWX",wxinit,1,con); wOp = shape(wxOp,2,ncol(wxOp)/2)[1,]; xOp = shape(wxOp,2,ncol(wxOp)/2)[2,]; wOp = (wOp##2)/sum(wOp##2); call Consolidate(xOp,wOp,1e-2,1e-5); if (ncol(xOp) ^= nrow(last)) then Improvement = 999; else Improvement = max(abs(last - (xOp` || wOp`))); end; /* / Table 18.2 /------------------------------------------------------------------*/ print "Table 18.2 - Continuous optimum designs",, "Region ------- Design ------- Crit",, "Convex [0,1] " ("x" // "w") ( (xOp || .) // (wOp || .)) [format=6.4] ' ' (exp(-DCrit(wOp,xOp))) [format=5.2],, "20-point Grid " ("x" // "w") ( (x20 || .) // (w20 || .)) [format=6.4] ' ' (exp(-DCrit(w20,x20))) [format=5.2],, "10-point Grid " ("x" // "w") ( x10 // w10 ) [format=6.4] ' ' (exp(-DCrit(w10,x10))) [format=5.2]; /* / Compute derivative function, equation (18.10) in the text. /------------------------------------------------------------------*/ start dd(xx,w,x) global(tPop,pPop); EMd = 0; EM = 0; free di; do i = 1 to ncol(tPop); fi = tPop[i]*( x#(1 - tPop[i]* x))#(0 < x)#( x < 1/tPop[i]); fx = tPop[i]*(xx*(1 - tPop[i]*xx))*(0 < xx)*(xx < 1/tPop[i]); Mi = fi*diag(w)*fi`; EMd = EMd + pPop[i]*det(inv(Mi))*fx`*inv(Mi)*fx; EM = EM + pPop[i]*det(inv(Mi)); end; return(EMd/EM); finish; free data; do xx = 0 to 1 by 0.01; data = data // (xx || dd(xx,w10,x10)); end; create fig1801 var ({x d}); append from data; title2 "Figure 18.1 - Derivative function"; %let pagesize = %sysfunc(getoption(pagesize)); options ps=40; proc plot data=fig1801; plot d*x; run; options ps=&pagesize; title2; title1;