*/**********************************************************************/ /* Atkinson, Donev, and Tobias, "Optimum Experimental Designs" */ /* */ /* NAME: Program 21.5 */ /* TITLE: Example 21.4: Two linear models */ /* KEYS: */ /* DATA: */ /* */ /* Author: Randy Tobias */ /* History: */ /* Created...............................................27Jun2007 */ /**********************************************************************/ title1 "Example 21.4: Two linear models"; data Can; do x = -1 to 1 by 0.05; x10 = 1; x11 = exp( x); x12 = exp(-x); x20 = 1; x21 = x; x22 = x*x; output; end; run; proc iml; use Can; read all var {x} into xCan; 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; x = shape(x,nrow(x)*ncol(x),1); w = shape(w,nrow(w)*ncol(w),1); w = w/sum(w); finish; start MakeZ(x); F1 = exp(x) || exp(-x); F2 = x || x#x ; return(j(nrow(x),1) || F1 || F2); finish; eps = 0; start DCritW(w) global(F,k,eps); p = ncol(F); r1 = 3; s1 = p - r1; r2 = 3; s2 = p - r2; F1 = F[,{1 2 3}]; F2 = F[,{1 4 5}]; Dw = diag((w##2)/sum(w##2)); m11 = F1`*Dw*F1 + eps*i(ncol(F1)); m22 = F2`*Dw*F2 + eps*i(ncol(F2)); M = F `*Dw*F + eps*i(ncol(F )); return( ((k/r1)*log(det(M11)) + ((1-k)/s1)*(log(det(M)) - log(det(M11)))) + ((k/r2)*log(det(M22)) + ((1-k)/s2)*(log(det(M)) - log(det(M22)))) ); finish; start DCritWX(wx) global(k,eps); n = nrow(wx)*ncol(wx)/2; w = shape(wx,2,n)[1,]`; w = (w##2)/sum(w##2); x = shape(wx,2,n)[2,]`; F1 = j(n,1) || exp(x) || exp(-x); F2 = j(n,1) || x || x#x; F = F1 || F2[,{2 3}]; p = ncol(F); r1 = 3; s1 = p - r1; r2 = 3; s2 = p - r2; Dw = diag((w##2)/sum(w##2)); m11 = F1`*Dw*F1 + eps*i(ncol(F1)); m22 = F2`*Dw*F2 + eps*i(ncol(F2)); M = F `*Dw*F + eps*i(ncol(F )); return( ((k/r1)*log(det(M11)) + ((1-k)/s1)*(log(det(M)) - log(det(M11)))) + ((k/r2)*log(det(M22)) + ((1-k)/s2)*(log(det(M)) - log(det(M22)))) ); finish; start DOptimize; x = XCan; F = MakeZ(x); winit = sqrt(j(nrow(X),1)/nrow(X)); call nlpqn(rc,wOp,"DCritW",winit,1); wOp = ((wOp##2)/sum(wOp##2))`; XOp = X [loc(wOp > 1e-6),]; wOp = wOp[loc(wOp > 1e-6),]; wOp = wOp /sum(wOp); x = xOp; F = MakeZ(x); winit = sqrt(wOp); call nlpqn(rc,wOp,"DCritW",winit,1); wOp = ((wOp##2)/sum(wOp##2))`; XOp = X [loc(wOp > 1e-6),]; wOp = wOp[loc(wOp > 1e-6),]; wOp = wOp /sum(wOp); eps = 1e-8; i = 0; diff = 1; do while ((i < 100) & (diff > 1e-8)); i = i + 1; last = (xOp || wOp); n = nrow(wOp)*ncol(wOp); con = (shape(.,1,n) || shape(-1,1,n)) // (shape(.,1,n) || shape( 1,1,n)); wxinit = sqrt(wOp)` || xOp`; call nlpqn(rc,wxOp,"DCritWX",wxinit,1,con); wOp = shape(wxOp,2,n)[1,]; wOp = (wOp##2)/sum(wOp##2); xOp = shape(wxOp,2,n)[2,]; call Consolidate(xOp,wOp,1e-5,1e-3); n = nrow(wOp)*ncol(wOp); if (n ^= nrow(last)) then diff = 1; else diff = max(abs(last - (xOp || wOp))); end; finish; k = 0.5; call DOptimize; wOp5 = wOp; xOp5 = xOp; k = 0.6; call DOptimize; wOp6 = wOp; xOp6 = xOp; start DCrit1WX(wx) global(eps); n = nrow(wx)*ncol(wx)/2; w = shape(wx,2,n)[1,]`; w = (w##2)/sum(w##2); x = shape(wx,2,n)[2,]`; F1 = j(n,1) || exp(x) || exp(-x); r1 = 3; Dw = diag((w##2)/sum(w##2)); m11 = F1`*Dw*F1 + eps*i(ncol(F1)); return(exp(log(det(M11))/r1)); finish; start DCrit2WX(wx) global(eps); n = nrow(wx)*ncol(wx)/2; w = shape(wx,2,n)[1,]`; w = (w##2)/sum(w##2); x = shape(wx,2,n)[2,]`; F2 = j(n,1) || x || x#x; r2 = 3; Dw = diag((w##2)/sum(w##2)); m22 = F2`*Dw*F2 + eps*i(ncol(F2)); return(exp(log(det(M22))/r2)); finish; %let t10 = 4.5; %let t11 = -1.5; %let t12 = -2.0; t1 = {&t10 &t11 &t12}`; start TCritW(_w) global(F1,F2,t1); w = shape(_w,nrow(_w)*ncol(_w),1); w = (w##2)/sum(w##2); Dw = diag(w); eta1 = F1*t1; eta2 = F2*inv(F2`*Dw*F2)*F2`*Dw*F1*t1; return(sum(w#((eta1 - eta2)##2))); finish; start TCritWX(wx) global(t1); n = nrow(wx)*ncol(wx)/2; w = shape(wx,2,n)[1,]`; w = (w##2)/sum(w##2); x = shape(wx,2,n)[2,]`; F1 = j(n,1) || exp(x) || exp(-x); F2 = j(n,1) || x || x#x; eta1 = F1*t1; eta2 = F2*inv(F2`*diag(w)*F2)*F2`*diag(w)*F1*t1; return(sum(w#((eta1 - eta2)##2))); finish; x = XCan; n = nrow(x); F1 = j(n,1) || exp(x) || exp(-x); F2 = j(n,1) || x || x#x; winit = sqrt(j(nrow(X),1)/nrow(X)); call nlpqn(rc,wOp,"TCritW",winit,1); xOp = x; call Consolidate(xOp,wOp,0,1e-4); eps = 1e-8; i = 0; diff = 1; do while ((i < 100) & (diff > 1e-8)); i = i + 1; last = (xOp || wOp); n = nrow(wOp)*ncol(wOp); con = (shape(.,1,n) || shape(-1,1,n)) // (shape(.,1,n) || shape( 1,1,n)); wxinit = sqrt(wOp)` || xOp`; call nlpqn(rc,wxOp,"TCritWX",wxinit,1,con); wOp = shape(wxOp,2,n)[1,]; wOp = (wOp##2)/sum(wOp##2); xOp = shape(wxOp,2,n)[2,]; call Consolidate(xOp,wOp,1e-2,1e-2); n = nrow(wOp)*ncol(wOp); if (n ^= nrow(last)) then diff = 1; else diff = max(abs(last - (xOp || wOp))); end; xOpT = xOp; wOpT = wOp; wOp12 = { 1 1 1}`/3; xOp12 = {-1 0 1}`; wxOp12 = wOp12 // xOp12; wxOpT = wOpT // xOpT ; wxOp5 = wOp5 // xOp5 ; wxOp6 = wOp6 // xOp6 ; Support = putn((xOp5+xOp6)`/2,'6.3'); print "Table 21.5 - Individual model and T-efficiencies for designs",, (wOp5`) [format=6.3 colname=Support rowname="w=0.5"] " " (DCrit1WX(wxOp5)/DCrit1WX(wxOp12)) [format=5.3 colname="EffD1"] (DCrit2WX(wxOp5)/DCrit2WX(wxOp12)) [format=5.3 colname="EffD2"] ( TCritWX(wxOp5)/ TCritWX(wxOpT )) [format=5.3 colname="EffT" ], (wOp6`) [format=6.3 rowname="w=0.6"] " " (DCrit1WX(wxOp6)/DCrit1WX(wxOp12)) [format=5.3] (DCrit2WX(wxOp6)/DCrit2WX(wxOp12)) [format=5.3] ( TCritWX(wxOp6)/ TCritWX(wxOpT )) [format=5.3]; title1;