options double crt memory=8;
name ANOVA 'NIST StRD ANOVA (via OLS) benchmarks';
const it,0 nmodel,11;
load;
mform(nrow=nmodel,ncol=1) tab;  ? Fstat
supres smpl nob coef;

? Si Resistivity
const nt,5 nrep,5;
const cfstat 1.18046237440255E+00 ;
set nob = nt*nrep;
smpl 1,nob;
read y;
196.3052 196.1240 196.1890 196.2569 196.3403 196.3042 196.3825 196.1669
196.3257 196.0422 196.1303 196.2005 196.2889 196.0343 196.1811 196.2795
196.1748 196.1494 196.1485 195.9885 196.2119 196.1051 196.1850 196.0052
196.2090
cert;

? Simon-Lesage 1,2,3
const nt,9 nrmax,2001;
set nob = nt*nrmax;
smpl 1,nob;
ifrac = 0;
trend(per=nrmax) obsi;
? 1st treatment
smpl 1,2001;
iodd = 2*int(obsi/2) .ne. obsi;  ? dummy for odd observations
ifrac = 3 + 2*iodd;  ? 5 3 5 3 5 ... 3 5
set ifrac(1) = 4;    ? 4 3 5 3 5 ... 3 5
? 2nd treatment  (also 4 6 8)
smpl 2002,4002;
ifrac = ifrac(-2001) - 1;  ? 3 2 4 2 4 ... 2 4
? 3rd treatment  (also 5 7 9)
smpl 4003,6003;
ifrac = ifrac(-2001) + 2;  ? 5 4 6 4 6 ... 4 6
? treatments 4-9 (same as 2-3)
smpl 6004,nob;
ifrac = ifrac(-4002);
smpl 1,nob;
y = 1 + ifrac/10;  ? Simon-Lesage 1,2,3
const nrep,21;   ? Simon-Lesage 1
set cfstat=nrep;
select obsi <= nrep;
cert;
const nrep,201;   ? Simon-Lesage 2
set cfstat=nrep;
select obsi <= nrep;
cert;
const nrep,2001;   ? Simon-Lesage 3
set cfstat=nrep;
select obsi <= nrep;
cert;

? Atomic Weight of Silver
const nt,2 nrep,24;
const cfstat 1.59467335677930E+01 ;
set nob = nt*nrep;
smpl 1,nob;
read y;
107.8681568 107.8681465 107.8681572 107.8681785 107.8681446 107.8681903
107.8681526 107.8681494 107.8681616 107.8681587 107.8681519 107.8681486
107.8681419 107.8681569 107.8681508 107.8681672 107.8681385 107.8681518
107.8681662 107.8681424 107.8681360 107.8681333 107.8681610 107.8681477
107.8681079 107.8681344 107.8681513 107.8681197 107.8681604 107.8681385
107.8681642 107.8681365 107.8681151 107.8681082 107.8681517 107.8681448
107.8681198 107.8681482 107.8681334 107.8681609 107.8681101 107.8681512
107.8681469 107.8681360 107.8681254 107.8681261 107.8681450 107.8681368
cert;

? Simon-Lesage 4-6
const nt,9 nrmax,2001;
set nob = nt*nrmax;
smpl 1,nob;
y = 1000000 + ifrac/10;  ? Simon-Lesage 4-6
const nrep,21;   ? Simon-Lesage 4
set cfstat=nrep;
select obsi <= nrep;
cert;
const nrep,201;   ? Simon-Lesage 5
set cfstat=nrep;
select obsi <= nrep;
cert;
const nrep,2001;   ? Simon-Lesage 6
set cfstat=nrep;
select obsi <= nrep;
cert;

? Simon-Lesage 7-9
const nt,9 nrmax,2001;
set nob = nt*nrmax;
smpl 1,nob;
y = 1000000000000 + ifrac/10;  ? Simon-Lesage 7-9
const nrep,21;   ? Simon-Lesage 7
set cfstat=nrep;
select obsi <= nrep;
cert;
const nrep,201;   ? Simon-Lesage 8
set cfstat=nrep;
select obsi <= nrep;
cert;
const nrep,2001;   ? Simon-Lesage 9
set cfstat=nrep;
select obsi <= nrep;
cert;

write(format="(/' Anova results'/
                9X,' NDigits'/
       ' Model   ','  Fstat  Pass=1/Fail=0')");
smpl 1,nmodel;
unmake tab c1;
c2 = (c1 >= 7);  ? Pass/Fail criteria
mat npass = sum(c2);
mmake tab cnames c1-c2;
write(format="(1X,A8,F7.1,F5.0)") tab;
write(format="(' Anova summary:  ',F3.0,'/11 passed.')") npass;

?======================================================
Proc Cert;
? Procedure for running NIST StRD Anova benchmarks --
? Record number of correct digits in table.
? Inputs (global variables)
?  nt   = number of treatments (number of dummies)
?  nrep = number of replicates (number of obs per cell)
? cfstat = certified F-statistic
?  tab  = table of output results
?  it   = pointer to row of tab for this model (maintained here)
?  by Clint Cummins  9/97

  set it = it+1;   ? pointer to row of tables
?
? create the dummy variables
?
  trend obs;
  treat = 1 + int( (obs-1)/nrep );  ? = 1 for first nrep obs, then 2, etc.
  list(prefix=dum,last=nt) dums;
  dummy treat dums;
  olsq(silent) y dums;
  ndigits @fst cfstat nd;  set tab(it,1) = nd;
endproc;
?======================================================
Proc ndigits try true nd;
? Count the number of digits to which 2 numbers agree
  local dc;
  set dc = abs[(try-true)/
               (true + (true=0))];  ? denom is replaced by 1 if true=0
  set dc = (dc >= 1e-15)*dc + (dc < 1e-15)*1e-15;  ? avoid log of zero
  set nd = -log10(dc);
endproc;
?======================================================
end;
noprint;
smpl 1,nmodel;
read(format='(A8)') names;
SiResist
Sim-Les1
Sim-Les2
Sim-Les3
AgAtomic
Sim-Les4
Sim-Les5
Sim-Les6
Sim-Les7
Sim-Les8
Sim-Les9
mmake cnames names;