/**********************************************************************/ /* Atkinson, Donev, and Tobias, "Optimum Experimental Designs" */ /* */ /* NAME: Program 17.5 */ /* TITLE: Example 17.5 - Sequential nonlinear optimum design */ /* KEYS: */ /* DATA: */ /* */ /* Author: Randy Tobias */ /* History: */ /* Created...............................................27Jun2007 */ /**********************************************************************/ title1 "Example 17.5 - Sequential nonlinear optimum design"; title2 "True optimum design"; data Can; do x1 = 0 to 3 by .05; do x2 = 0 to 3 by .05; output; end; end; run; %let Truet1 = 2.9; %let Truet2 = 12.2; %let Truet3 = 0.69; data Jac; set Can; d1 = &Truet3*x1*(1 + &Truet2*x2) / ((1 + &Truet1*x1 + &Truet2*x2)**2); d2 = -&Truet1*&Truet3*x1*x2 / ((1 + &Truet1*x1 + &Truet2*x2)**2); d3 = &Truet1*x1 / (1 + &Truet1*x1 + &Truet2*x2); run; proc optex data=Jac coding=none noprint; model d1 d2 d3 / noint; generate criterion=d n=3 method=m_fedorov niter=1000 keep=10; output out=OptDesign; id x1 x2; proc print data=OptDesign noobs; var x1 x2; run; title2; /* / A macro that implements the iterative procedure described in / Section 17.7 of the text. Starting with an initial design, we / iterate three steps: / * Optimal augmentation / * Simulation / * Fit and check / until the maximum relative width of the parameter confidence / is intervals less than LIMIT. /---------------------------------------------------------------------*/ %macro Add(Limit=0.1,nAdd=1,TrueSig=0.01,seed=1); %local notesopt; %let notesopt = %sysfunc(getoption(notes)); options nonotes; /* / The assumed true parameters /------------------------------------------------------------------*/ %let Truet1 = 2.9; %let Truet2 = 12.2; %let Truet3 = 0.69; /* / Initialize parameter estimates by fitting model to initial / points. /------------------------------------------------------------------*/ ods listing close; proc nlin data=InitDesign; model y = t1*t3*x1 / (1 + t1*x1 + t2*x2); parms t1=10.39 t2=48.83 t3=0.74; ods output ParameterEstimates=Parm; run; ods listing; data _null_; set Parm; call symput('Fit'||trim(left(Parameter)),trim(left(Estimate))); run; data Experiment; set InitDesign; ey = .; e = .; nNow = 0; run; data _null_; set Experiment; call symput('nD',_N_); run; /* / Fit and check step: Fit current data, update parameters, and / check the relative CL widths to see if we're done. /------------------------------------------------------------------*/ ods listing close; proc nlin data=Experiment outest=Parm; parm t1=&Fitt1 t2=&Fitt2 t3=&Fitt3; model y = t1*t3*x1 / (1 + t1*x1 + t2*x2); ods output ParameterEstimates=ParmCL; run; ods listing; data _null_; set Parm(where=(_TYPE_='FINAL')); call symput('Fitt1',trim(left(t1))); call symput('Fitt2',trim(left(t2))); call symput('Fitt3',trim(left(t3))); data ParmCL; set ParmCL; RelWidth = (UpperCL - LowerCL) / abs(Estimate); proc summary data=ParmCL; var RelWidth; output out=maxRelWidth max=maxRelWidth; data _null_; set maxRelWidth; call symput('Done',trim(left(maxRelWidth < &Limit))); run; %let nLast = &nD; %do %while (^&Done); /* / Optimal augmentation step: Compute the Jacobian around the / current parameters, for both the current design and for the / candidates, and optimally augment the former with nAdd points / from the latter. /---------------------------------------------------------------*/ data LastDesign; set Experiment; d1 = &Fitt3*x1*(1 + &Fitt2*x2) / ((1 + &Fitt1*x1 + &Fitt2*x2)**2); d2 = -&Fitt1*&Fitt3*x1*x2 / ((1 + &Fitt1*x1 + &Fitt2*x2)**2); d3 = &Fitt1*x1 / (1 + &Fitt1*x1 + &Fitt2*x2); data Jac; set Can; d1 = &Fitt3*x1*(1 + &Fitt2*x2) / ((1 + &Fitt1*x1 + &Fitt2*x2)**2); d2 = -&Fitt1*&Fitt3*x1*x2 / ((1 + &Fitt1*x1 + &Fitt2*x2)**2); d3 = &Fitt1*x1 / (1 + &Fitt1*x1 + &Fitt2*x2); run; %let nD = %eval(&nD + &nAdd); data LastDesign; set LastDesign; data Jac; set Jac; y = .; ey = .; e = .; nNow = .; proc optex data=Jac coding=none noprint; model d1 d2 d3 / noint; generate criterion=d n=&nD method=m_fedorov augment=LastDesign; output out=Design; id x1 x2 y ey e nNow; data _null_; set Design; call symput('nNow',trim(left(_N_))); run; /* / Simulation step: Generate response values for the new points. /---------------------------------------------------------------*/ data Experiment; set Design; if (y = .) then do; ey = &Truet1*&Truet3*x1 / (1 + &Truet1*x1 + &Truet2*x2); e = &TrueSig*rannor(&seed); call symput('seed',trim(left((2**31-1)*ranuni(&seed)))); y = ey + e; nNow = &nNow; end; run; /* / Fit and check step: Fit current data, update parameters, and / check the relative CL widths to see if we're done. /---------------------------------------------------------------*/ ods listing close; proc nlin data=Experiment outest=Parm; parm t1=&Fitt1 t2=&Fitt2 t3=&Fitt3; model y = t1*t3*x1 / (1 + t1*x1 + t2*x2); ods output ParameterEstimates=ParmCL; run; ods listing; data _null_; set Parm(where=(_TYPE_='FINAL')); call symput('Fitt1',trim(left(t1))); call symput('Fitt2',trim(left(t2))); call symput('Fitt3',trim(left(t3))); data ParmCL; set ParmCL; RelWidth = (UpperCL - LowerCL) / abs(Estimate); proc summary data=ParmCL; var RelWidth; output out=maxRelWidth max=maxRelWidth; data _null_; set maxRelWidth; call symput('Done',trim(left(maxRelWidth < &Limit))); run; %end; options ¬esopt; %mend; data InitDesign; input x1 x2 y; cards; 1 1 0.126 2 1 0.219 1 2 0.076 2 2 0.126 ; %Add(seed=76); proc sort data=Experiment; by nNow x2 x1; run; title2 "Sequential design, converging to true optimum"; proc print data=Experiment; var nNow x1 x2 y; run; title2; title1;