/**********************************************************************/ /* Atkinson, Donev, and Tobias, "Optimum Experimental Designs" */ /* */ /* NAME: Program 11.1 */ /* TITLE: Example 11.1 - Cubic regression through the origin */ /* KEYS: */ /* DATA: */ /* */ /* Author: Randy Tobias */ /* History: */ /* Created...............................................27Jun2007 */ /**********************************************************************/ title1 "Table 11.1 & Figure 11.5 - Cubic regression through the origin"; title2 "Sequential optimization"; /* / There are two ways to build the 50-point design by sequentially / adding points that maximize d(x). The SEQUENTIAL method in / OPTEX does precisely this, but we must remember to turn off the / automatic coding: optimal designs for typical linear models with / intercepts are invariant to coding, but this design is not. /---------------------------------------------------------------------*/ data InitDesign; input x @@; output; datalines; 0.1 0.5 0.9 ; data Candidates; do x = 0 to 1 by 0.05; output; end; run; proc optex data=Candidates coding=none; model x x*x x*x*x / noint; generate n=50 augment=InitDesign method=sequential; output out=SeqDesign; run; /* / The other way is to build the sequential design "by hand", as it / were. This way also allows us to see the first four points added. /---------------------------------------------------------------------*/ proc iml; /* / Define the prediction grid. /---------------------------------------------------------------------*/ use Candidates; read all var {"x"} into xGrid; Z = xGrid || xGrid#xGrid || xGrid#xGrid#xGrid; /* / The standardized variance function. /---------------------------------------------------------------------*/ start d(x) global(Z); F = x || x#x || x#x#x; V = inv(F`*F); free d; do i = 1 to nrow(Z); d = d // Z[i,]*V*Z[i,]`*nrow(x); end; return(d); finish; free Additions; /* / Beginning with an arbitrary design ... /---------------------------------------------------------------------*/ use InitDesign; read all var {"x"}; /* / ... compute the standardized variance over the prediction grid, ... /---------------------------------------------------------------------*/ d = d(x); /* / ... and determine the point which maximizes the standardized / variance. /---------------------------------------------------------------------*/ xAdd = xGrid[d[<:>],1]; Additions = Additions // (nrow(x)-3 || xAdd || d[<>]); /* / Add the maximizing point to the design and repeat. /---------------------------------------------------------------------*/ x = x // xAdd; d = d(x); xAdd = xGrid[d[<:>]]; Additions = Additions // (nrow(x)-3 || xAdd || d[<>]); x = x // xAdd; d = d(x); xAdd = xGrid[d[<:>]]; Additions = Additions // (nrow(x)-3 || xAdd || d[<>]); x = x // xAdd; d = d(x); xAdd = xGrid[d[<:>]]; Additions = Additions // (nrow(x)-3 || xAdd || d[<>]); print "The first four additions",,Additions; do while (nrow(x) < 50); x = x // xAdd; d = d(x); xAdd = xGrid[d[<:>]]; Additions = Additions // (nrow(x)-3 || xAdd || d[<>]); end; create Design var {"x"}; append from x; quit; title3 "The first 50 points"; proc freq data=Design; table x / nocum; run; title3; title3 "The two sequential methods match"; proc sort data=SeqDesign; by x; proc sort data= Design; by x; proc compare data=SeqDesign compare=Design brief; run; title3; title2; /* / The continuous optimization proceeds very much as in Program 9.1, / except that the model has changed and the D-criterion is being / optimized. Notice that the values of the D-criterion |M| itself / are typically too small for a numerical optimizer to detect / changes, so we actually use log(|M|) instead. /---------------------------------------------------------------------*/ title2 "Continuous optimization"; options nonotes; proc iml; start Crit(w,x); Z = x || x#x || x#x#x; Mi = Z`*diag(w)*Z; return(log(det(Mi))); finish; start CritW(w) global(x); return(Crit((w##2)/sum(w##2),x)); finish; start CritX(x) global(w); return(Crit(w,x)); finish; start CritWX(wx); w = shape(wx,2,ncol(wx)/2)[1,]`; w = (w##2)/sum(w##2); x = shape(wx,2,ncol(wx)/2)[2,]`; return(Crit(w,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 (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); x = xNew; w = wNew; x[r] = xNew; w[r] = wNew; finish; start tlc(n); return(trim(left(char(n,12)))); finish; %let xLo = 0; %let xHi = 1; %let Inc = 40; x = &xLo + (&xHi - &xLo)*(0:&Inc)`/&Inc; winit = j(nrow(x),1) / nrow(x); call nlpqn(rc,wOp,"CritW",winit,1); wOp = (wOp##2)/sum(wOp##2); xOp = x`; print "Initial optimum weights over support grid",, (xOp`) [colname={"Support"} format=7.4] (wOp`) [colname={"Weights"} format=7.4] (Crit(wOp`,xOp`)) [colname={"log(|M|)"} format=6.4]; call Consolidate(xOp,wOp,0,1e-4); print "Simplified inital optimum design",, (xOp`) [colname={"Support"} format=7.4] (wOp`) [colname={"Weights"} format=7.4] (Crit(wOp`,xOp`)) [colname={"log(|M|)"} format=6.4]; 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(&xLo,1,ncol(xOp))) // (shape(.,1,ncol(xOp)) || shape(&xHi,1,ncol(xOp))); call nlpqn(rc,wxOp,"CritWX",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`))); print ("Improving both weight and support: iteration " + tlc(i)),, (xOp`) [colname={"Support"} format=7.4] (wOp`) [colname={"Weights"} format=7.4] (Crit(wOp`,xOp`)) [colname={"log(|M|)/3"} format=6.4] Improvement; end; wOp = wOp`; xOp = xOp`; print "Final optimum design",, (1*xOp) [colname={"Support"} format=7.4] (1*wOp) [colname={"Weights"} format=7.4] (Crit(wOp,xOp)) [colname={"log(|M|)/3"} format=6.4]; /* / Finish by checking the maximum variance of prediction over the / design region, which should be the same as the number of / parameters, namely, 3. /---------------------------------------------------------------------*/ free xGrid; do x = 0 to 1 by 0.01; xGrid = xGrid // x; end; Z = xGrid || xGrid#xGrid || xGrid#xGrid#xGrid; start d(x,w) global(Z); F = x || x#x || x#x#x; V = inv(F`*diag(w)*F); free d; do i = 1 to nrow(Z); d = d // Z[i,]*V*Z[i,]`; end; return(d); finish; d = d(xOp,wOp); print "Maximum standardized variance",,(d[<>]); title2; title1;