/**********************************************************************/ /* Atkinson, Donev, and Tobias, "Optimum Experimental Designs" */ /* */ /* NAME: Program 18.3 */ /* TITLE: Tables 18.5 - 18.7: Exponential decay, optimum designs */ /* for various priors */ /* KEYS: */ /* DATA: */ /* NOTES: Took ~40 seconds to run on server-class PC in May, 2007. */ /* */ /* Author: Randy Tobias */ /* History: */ /* Created...............................................27Jun2007 */ /**********************************************************************/ title1 "Example 18.2 - Exponential decay: optimum designs for various priors"; proc iml; /* / 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; start llog(x); if (x > 0) then return(log(x)); else return(-9999); finish; 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 = x#exp(-tPop[i]*x); 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; /* / As Figure 18.2 in the text indicates, this optimization problem is / very flat, with many local optima. To improve the likelihood of / finding the true optimum, we use many random initial designs. /---------------------------------------------------------------------*/ start Optimize(xinit) global(x); niter = 20; ninit = ncol(xinit); lLo = log(min(xinit)); lHi = log(max(xinit)); do iter = 0 to niter; if (iter = 0) then x = xinit; else x = exp(lLo + (lHi - lLo)*ranuni(1:ninit)); winit = j(ninit,1) / ninit; call nlpqn(rc,wOp,"DCritW",winit,1); wOp = ((wOp##2)/sum(wOp##2)); xOp = x; call Consolidate(xOp,wOp,0,1e-5); i = 0; Improvement = 1; do while ((i < 100) & (Improvement > 1e-6)); i = i + 1; last = (xOp` || wOp`); wxinit = sqrt(wOp) || xOp; con = (shape(.,1,ncol(xOp)) || shape( 0,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-3,1e-5); if (ncol(xOp) ^= nrow(last)) then Improvement = 999; else Improvement = max(abs(last - (xOp` || wOp`))); end; dOp = DCrit(wOp,xOp); if (iter = 0) then do; wBest = wOp; xBest = xOp; dBest = dOp; end; else if (dOp > dBest) then do; wBest = wOp; xBest = xOp; dBest = dOp; end; end; return(wBest // xBest); finish; options nonotes; Criterion = 1; /* / Find optimum criterion I designs for different log-uniform priors. /---------------------------------------------------------------------*/ nus = {1 3 7 13 100}; do inu = 1 to ncol(nus); nu = nus[inu]; tPop = 1/nu || 1/sqrt(nu) || 1 || sqrt(nu) || nu; pPop = j(1,ncol(tPop))/ncol(tPop); nint = 20; wxOp = Optimize(max(tPop)*2**(-(0:nint))); wOp = wxOp[1,]; xOp = wxOp[2,]; if (inu=1) then print "Table 18.5 - Dependence of design on range of prior distribution"; Support = xOp`; Weights = wOp`; print nu Support [format=6.4] Weights [format=6.4]; if (nu = 1) then do; xOp1 = xOp; wOp1 = wOp; end; if (nu = 7) then do; xOp2 = xOp; wOp2 = wOp; xOp3 = xOp; wOp3 = {2 1 1}/4; end; end; /* / Find optimum criterion I designs for log-normal priors with / different numbers of points. The designs aren't exactly the same / as in the text, so we also print the critical value to compare. /---------------------------------------------------------------------*/ tPop7 = 1/7 || 1/sqrt(7) || 1 || sqrt(7) || 7; v7 = sqrt(ssq(log(tPop7) - (log(tPop7))[:])/ncol(tPop7)); npts = {7 15}; do inpt = 1 to ncol(npts); npt = npts[inpt]; lnPop = -2.5 + (0:npt)*(2.5 - -2.5)/npt; free tPop pPop; do i = 1 to ncol(lnPop)-1; tPop = tPop || exp(( v7*lnPop[i+1] + v7*lnPop[i] )/2); pPop = pPop || probnorm( lnPop[i+1]) - probnorm( lnPop[i]); end; pPop = pPop / sum(pPop); nint = 20; wxOp = Optimize(max(tPop)*2**(-(0:nint))); wOp = wxOp[1,]; xOp = wxOp[2,]; if (inpt=1) then print "Table 18.6 - Optimum designs for discretized lognormal priors"; if (npt = 7) then do; xADT = {0.1013 0.2295 0.6221 1.6535 4.2724}; wADT = {0.0947 0.1427 0.3623 0.2549 0.1454}; end; else do; xADT = {0.1077 0.3347 0.7480 1.4081 3.7769}; wADT = {0.1098 0.2515 0.2153 0.2491 0.1743}; end; Support = xOp`; Weights = wOp`; DCritOp = DCrit(wOp,xOp); DCritADT = DCrit(wADT,xADT); print npt Support [format=6.4] Weights [format=6.4] DCritOp DCritADT; if (npt = 7) then do; xOp4 = xOp; wOp4 = wOp; end; if (npt = 15) then do; xOp5 = xOp; wOp5 = wOp; end; end; /* / Evaluate designs against the one that's optimum for the 15-point / log-normal prior. /---------------------------------------------------------------------*/ npt = 15; lnPop = -2.5 + (0:npt)*(2.5 - -2.5)/npt; free tPop pPop; do i = 1 to ncol(lnPop)-1; tPop = tPop || exp(( v7*lnPop[i+1] + v7*lnPop[i] )/2); pPop = pPop || probnorm( lnPop[i+1]) - probnorm( lnPop[i]); end; pPop = pPop / sum(pPop); dOp1 = exp(DCrit(wOp1,xOp1)); dOp2 = exp(DCrit(wOp2,xOp2)); dOp3 = exp(DCrit(wOp3,xOp3)); dOp4 = exp(DCrit(wOp4,xOp4)); dOp5 = exp(DCrit(wOp5,xOp5)); Design = ( "1 point " // "5 point log-uniform, nu = 7" // "Exact N=4, nu=7 " // "7 point log-normal " // "15 point log-normal "); Efficiency = ((dOp1 // dOp2 // dOp3 // dOp4 // dOp5)/dOp5); print "Table 18.7 - Efficiencies for optimum designs for various priors",, Design Efficiency [format=percent10.2]; title1;