#delim ;
*
 Certify powercal using generated data.
*;

* Begin certification script *;
cscript "Certification of powercal"
 adofiles powercal;
version;

*
 Create data set with wide ranges of original values
 for the 5 parameters power, alpha, delta, sdinf and nunit,
 stored in 5 variables
 named powerold, alphaold, deltaold, sdinfold and nunitold,
 and set local macros for mid-values for each of these parameters.
 This allows the user to plot each ordered pair of parameters,
 with original values for the first parameter in the pair
 and values calculated by powercal for the second parameter in the pair,
 assuming mid-values for the other 3 parameters.
 The user can also check that the formula
 deriving the second parameter from the first
 is approximately the inverse of the formula
 deriving the first parameter from the second.
*;
local nobsmax=2001;
set obs `nobsmax';
gene int seqnum=_n;
* Mid-values for power, alpha, delta, sdinf and nunit *;
local powermid=0.5;
local alphamid=0.05;
local deltamid=2;
local sdinfmid=10;
local nunitmid=100;
* Other values to be tested for power, alpha, delta, sdinf and nunit *;
* Power *;
gene powerold=(_n-1)/(_N-1);
* Alpha *;
local alphamin=0.5e-6;
gene alphaold=exp( log(`alphamin')*(_n-1)/(_N-1) );
* Delta and Sdinf *;
local mmrat=8;
local deltamin=0;local deltamax=5*`deltamid';
gene deltaold=`deltamin' + (`deltamax'-`deltamin')*(_n-1)/(_N-1);
gene sdinfold=sqrt(`nunitmid')/deltaold;
* Nunit *;
gene int nunitold=_n-1;
* Variable labels *;
lab var seqnum "Sequence number";
foreach X in power alpha delta sdinf nunit {;lab var `X'old "Original `X'";};
* Describe data *;
desc;summ;


*
 Test sequence 1.
 Error conditions for powercal.
*;
* All input expressions provided *;
rcof "noisily powercal thingnew, p(`powermid') a(`alphamid') d(`deltamid') s(`sdinfmid') n(`nunitmid')"==498;


*
 Test sequence 2.
 Missing input values.
*;
* Calculation of power *;
powercal thingnew, a(.) d(`deltamid') s(`sdinfmid') n(`nunitmid');
assert missing(thingnew);drop thingnew;
powercal thingnew, a(`alphamid') d(.) s(`sdinfmid') n(`nunitmid');
assert missing(thingnew);drop thingnew;
powercal thingnew, a(`alphamid') d(`deltamid') s(.) n(`nunitmid');
assert missing(thingnew);drop thingnew;
powercal thingnew, a(`alphamid') d(`deltamid') s(`sdinfmid') n(.);
assert missing(thingnew);drop thingnew;
* Calculation of alpha *;
powercal thingnew, p(.) d(`deltamid') s(`sdinfmid') n(`nunitmid');
assert missing(thingnew);drop thingnew;
powercal thingnew, p(`powermid') d(.) s(`sdinfmid') n(`nunitmid');
assert missing(thingnew);drop thingnew;
powercal thingnew, p(`powermid') d(`deltamid') s(.) n(`nunitmid');
assert missing(thingnew);drop thingnew;
powercal thingnew, p(`powermid') d(`deltamid') s(`sdinfmid') n(.);
assert missing(thingnew);drop thingnew;
* Calculation of delta *;
powercal thingnew, p(.) a(`alphamid') s(`sdinfmid') n(`nunitmid');
assert missing(thingnew);drop thingnew;
powercal thingnew, p(`powermid') a(.) s(`sdinfmid') n(`nunitmid');
assert missing(thingnew);drop thingnew;
powercal thingnew, p(`powermid') a(`alphamid') s(.) n(`nunitmid');
assert missing(thingnew);drop thingnew;
powercal thingnew, p(`powermid') a(`alphamid') s(`sdinfmid') n(.);
assert missing(thingnew);drop thingnew;
* Calculation of sdinf *;
powercal thingnew, p(.) a(`alphamid') d(`deltamid') n(`nunitmid');
assert missing(thingnew);drop thingnew;
powercal thingnew, p(`powermid') a(.) d(`deltamid') n(`nunitmid');
assert missing(thingnew);drop thingnew;
powercal thingnew, p(`powermid') a(`alphamid') d(.) n(`nunitmid');
assert missing(thingnew);drop thingnew;
powercal thingnew, p(`powermid') a(`alphamid') d(`deltamid') n(.);
assert missing(thingnew);drop thingnew;
* Calculation of nunit *;
powercal thingnew, p(.) a(`alphamid') d(`deltamid') s(`sdinfmid');
assert missing(thingnew);drop thingnew;
powercal thingnew, p(`powermid') a(.) d(`deltamid') s(`sdinfmid');
assert missing(thingnew);drop thingnew;
powercal thingnew, p(`powermid') a(`alphamid') d(.) s(`sdinfmid');
assert missing(thingnew);drop thingnew;
powercal thingnew, p(`powermid') a(`alphamid') d(`deltamid') s(.);
assert missing(thingnew);drop thingnew;


*
 Test sequence 3.
 Out-of-range input values.
*;
* Add a variable containing a list of possibly out of range test inputs *;
gene testinput=_n-int(_N/2)+1;
lab var testinput "Test input";
* Calculation of power *;
powercal thingnew, a(testinput) d(`deltamid') s(`sdinfmid') n(`nunitmid');
assert !missing(thingnew)==((testinput>0)&(testinput<1));
drop thingnew;
powercal thingnew, a(`alphamid') d(testinput) s(`sdinfmid') n(`nunitmid');
assert !missing(thingnew)==((testinput>0)&(!missing(testinput)));
drop thingnew;
powercal thingnew, a(`alphamid') d(`deltamid') s(testinput) n(`nunitmid');
assert !missing(thingnew)==((testinput>0)&(!missing(testinput)));
drop thingnew;
powercal thingnew, a(`alphamid') d(`deltamid') s(`sdinfmid') n(testinput);
assert !missing(thingnew)==((testinput>0)&(!missing(testinput)));
drop thingnew;
* Calculation of alpha *;
powercal thingnew, p(testinput) d(`deltamid') s(`sdinfmid') n(`nunitmid');
assert !missing(thingnew)==((testinput>0)&(testinput<1));
drop thingnew;
powercal thingnew, p(`powermid') d(testinput) s(`sdinfmid') n(`nunitmid');
assert !missing(thingnew)==((testinput>0)&(!missing(testinput)));
drop thingnew;
powercal thingnew, p(`powermid') d(`deltamid') s(testinput) n(`nunitmid');
assert !missing(thingnew)==((testinput>0)&(!missing(testinput)));
drop thingnew;
powercal thingnew, p(`powermid') d(`deltamid') s(`sdinfmid') n(testinput);
assert !missing(thingnew)==((testinput>0)&(!missing(testinput)));
drop thingnew;
* Calculation of delta *;
powercal thingnew, p(testinput) a(`alphamid') s(`sdinfmid') n(`nunitmid');
assert !missing(thingnew)==((testinput>0)&(testinput<1));
drop thingnew;
powercal thingnew, p(`powermid') a(testinput) s(`sdinfmid') n(`nunitmid');
assert !missing(thingnew)==((testinput>0)&(testinput<1));
drop thingnew;
powercal thingnew, p(`powermid') a(`alphamid') s(testinput) n(`nunitmid');
assert !missing(thingnew)==((testinput>0)&(!missing(testinput)));
drop thingnew;
powercal thingnew, p(`powermid') a(`alphamid') s(`sdinfmid') n(testinput);
assert !missing(thingnew)==((testinput>0)&(!missing(testinput)));
drop thingnew;
* Calculation of sdinf *;
powercal thingnew, p(testinput) a(`alphamid') d(`deltamid') n(`nunitmid');
assert !missing(thingnew)==((testinput>0)&(testinput<1));
drop thingnew;
powercal thingnew, p(`powermid') a(testinput) d(`deltamid') n(`nunitmid');
assert !missing(thingnew)==((testinput>0)&(testinput<1));
drop thingnew;
powercal thingnew, p(`powermid') a(`alphamid') d(testinput) n(`nunitmid');
assert !missing(thingnew)==((testinput>0)&(!missing(testinput)));
drop thingnew;
powercal thingnew, p(`powermid') a(`alphamid') d(`deltamid') n(testinput);
assert !missing(thingnew)==((testinput>0)&(!missing(testinput)));
drop thingnew;
* Calculation of nunit *;
powercal thingnew, p(testinput) a(`alphamid') d(`deltamid') s(`sdinfmid');
assert !missing(thingnew)==((testinput>0)&(testinput<1));
drop thingnew;
powercal thingnew, p(`powermid') a(testinput) d(`deltamid') s(`sdinfmid');
assert !missing(thingnew)==((testinput>0)&(testinput<1));
drop thingnew;
powercal thingnew, p(`powermid') a(`alphamid') d(testinput) s(`sdinfmid');
assert !missing(thingnew)==((testinput>0)&(!missing(testinput)));
drop thingnew;
powercal thingnew, p(`powermid') a(`alphamid') d(`deltamid') s(testinput);
assert !missing(thingnew)==((testinput>0)&(!missing(testinput)));
drop thingnew;
* Calculation of nunit with the -noceil- option *;
powercal thingnew, p(testinput) a(`alphamid') d(`deltamid') s(`sdinfmid') noceil;
assert !missing(thingnew)==((testinput>0)&(testinput<1));
drop thingnew;
powercal thingnew, p(`powermid') a(testinput) d(`deltamid') s(`sdinfmid') noceil;
assert !missing(thingnew)==((testinput>0)&(testinput<1));
drop thingnew;
powercal thingnew, p(`powermid') a(`alphamid') d(testinput) s(`sdinfmid') noceil;
assert !missing(thingnew)==((testinput>0)&(!missing(testinput)));
drop thingnew;
powercal thingnew, p(`powermid') a(`alphamid') d(`deltamid') s(testinput) noceil;
assert !missing(thingnew)==((testinput>0)&(!missing(testinput)));
drop thingnew;
* Drop test input variable *;
drop testinput;


*
 Test sequence 4.
 Test each pair of powercal expression options.
*;

*
 Power and Alpha
*;
* Normal distribution, alpha to power (and back) *;
powercal powernew, a(alphaold) d(`deltamid') s(`sdinfmid') n(`nunitmid');
powercal alphanew, p(powernew) d(`deltamid') s(`sdinfmid') n(`nunitmid');
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>=0 if !missing(`Z');assert `Z'<=1 if !missing(`Z');
};
assert reldif(alphanew,alphaold)<=1e-9 if !missing(alphanew);
version 7:graph powernew alphaold, c(L) ylab xlab yline(`powermid') xline(`alphamid') xlog xreverse
 t1("Calculations assuming Normal distribution");
more;
drop *new;
* Normal distribution, power to alpha (and back) *;
powercal alphanew, p(powerold) d(`deltamid') s(`sdinfmid') n(`nunitmid');
powercal powernew, a(alphanew) d(`deltamid') s(`sdinfmid') n(`nunitmid');
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>=0 if !missing(`Z');assert `Z'<=1 if !missing(`Z');
};
assert reldif(powernew,powerold)<=1e-9 if !missing(powernew);
version 7:graph powerold alphanew, c(L) ylab xlab yline(`powermid') xline(`alphamid') xlog xreverse
 t1("Calculations assuming Normal distribution");
more;
drop *new;
* Student t-distribution, alpha to power (and back) *;
powercal powernew, a(alphaold) d(`deltamid') s(`sdinfmid') n(`nunitmid') td(`nunitmid'-1);
powercal alphanew, p(powernew) d(`deltamid') s(`sdinfmid') n(`nunitmid') td(`nunitmid'-1);
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>=0 if !missing(`Z');assert `Z'<=1 if !missing(`Z');
};
assert reldif(alphanew,alphaold)<=1e-9 if !missing(alphanew);
version 7:graph powernew alphaold, c(L) ylab xlab yline(`powermid') xline(`alphamid') xlog xreverse
 t1("Calculations assuming Student's t-distribution");
more;
drop *new;
* Student t-distribution, power to alpha (and back) *;
powercal alphanew, p(powerold) d(`deltamid') s(`sdinfmid') n(`nunitmid') td(`nunitmid'-1);
powercal powernew, a(alphanew) d(`deltamid') s(`sdinfmid') n(`nunitmid') td(`nunitmid'-1);
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>=0 if !missing(`Z');assert `Z'<=1 if !missing(`Z');
};
assert reldif(powernew,powerold)<=1e-9 if !missing(powernew);
version 7:graph powerold alphanew, c(L) ylab xlab yline(`powermid') xline(`alphamid') xlog xreverse
 t1("Calculations assuming Student's t-distribution");
more;
drop *new;

*
 Power and Delta
*;
* Normal distribution, delta to power (and back) *;
powercal powernew, a(`alphamid') d(deltaold) s(`sdinfmid') n(`nunitmid');
powercal deltanew, p(powernew) a(`alphamid') s(`sdinfmid') n(`nunitmid');
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>0 if !missing(`Z');
};
assert powernew<=1 if !missing(powernew);
assert reldif(deltanew,deltaold)<1e-2 if !missing(deltanew);
version 7:graph powernew deltaold, c(L) ylab xlab yline(`powermid') xline(`deltamid')
 t1("Calculations assuming Normal distribution");
more;
drop *new;
* Normal distribution, power to delta (and back) *;
powercal deltanew, p(powerold) a(`alphamid') s(`sdinfmid') n(`nunitmid');
powercal powernew, a(`alphamid') d(deltanew) s(`sdinfmid') n(`nunitmid');
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>0 if !missing(`Z');
};
assert powernew<=1 if !missing(powernew);
assert reldif(powernew,powerold)<1e-9 if !missing(powernew);
version 7:graph powerold deltanew, c(L) ylab xlab yline(`powermid') xline(`deltamid')
 t1("Calculations assuming Normal distribution");
more;
drop *new;
* Student t-distribution, delta to power (and back) *;
powercal powernew, a(`alphamid') d(deltaold) s(`sdinfmid') n(`nunitmid') td(`nunitmid'-1);
powercal deltanew, p(powernew) a(`alphamid') s(`sdinfmid') n(`nunitmid') td(`nunitmid'-1);
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>0 if !missing(`Z');
};
assert powernew<=1 if !missing(powernew);
assert reldif(deltanew,deltaold)<1e-2 if !missing(deltanew);
version 7:graph powernew deltaold, c(L) ylab xlab yline(`powermid') xline(`deltamid')
 t1("Calculations assuming Student's t-distribution");
more;
drop *new;
* Student t-distribution, power to delta (and back) *;
powercal deltanew, p(powerold) a(`alphamid') s(`sdinfmid') n(`nunitmid') td(`nunitmid'-1);
powercal powernew, a(`alphamid') d(deltanew) s(`sdinfmid') n(`nunitmid') td(`nunitmid'-1);
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>0 if !missing(`Z');
};
assert powernew<=1 if !missing(powernew);
assert reldif(powernew,powerold)<=1e-9 if !missing(powernew);
version 7:graph powerold deltanew, c(L) ylab xlab yline(`powermid') xline(`deltamid')
 t1("Calculations assuming Student's t-distribution");
more;
drop *new;

*
 Power and Sdinf
*;
* Normal distribution, sdinf to power (and back) *;
powercal powernew, a(`alphamid') d(`deltamid') s(sdinfold) n(`nunitmid');
powercal sdinfnew, p(powernew) a(`alphamid') d(`deltamid') n(`nunitmid');
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>0 if !missing(`Z');
};
assert powernew<=1 if !missing(powernew);
assert reldif(sdinfnew,sdinfold)<9e-3 if !missing(sdinfnew);
version 7:graph powernew sdinfold, c(L) ylab xlab yline(`powermid') xline(`sdinfmid') xlog
 t1("Calculations assuming Normal distribution");
more;
drop *new;
* Normal distribution, power to sdinf (and back) *;
powercal sdinfnew, p(powerold) a(`alphamid') d(`deltamid') n(`nunitmid');
powercal powernew, a(`alphamid') d(`deltamid') s(sdinfnew) n(`nunitmid');
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>0 if !missing(`Z');
};
assert powernew<=1 if !missing(powernew);
assert reldif(powernew,powerold)<1e-9 if !missing(powernew);
* Note that sdinfnew has some extreme values *;
version 7:graph powerold sdinfnew, c(L) xlab yline(`powermid') xline(`sdinfmid') xlog
 t1("Calculations assuming Normal distribution");
more;
drop *new;
* Student t-distribution, sdinf to power (and back) *;
powercal powernew, a(`alphamid') d(`deltamid') s(sdinfold) n(`nunitmid') td(`nunitmid'-1);
powercal sdinfnew, p(powernew) a(`alphamid') d(`deltamid') n(`nunitmid') td(`nunitmid'-1);
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>0 if !missing(`Z');
};
assert powernew<=1 if !missing(powernew);
assert reldif(sdinfnew,sdinfold)<9e-3 if !missing(sdinfnew);
version 7:graph powernew sdinfold, c(L) ylab xlab yline(`powermid') xline(`sdinfmid') xlog
 t1("Calculations assuming Student's t-distribution");
more;
drop *new;
* Student t-distribution, power to sdinf (and back) *;
powercal sdinfnew, p(powerold) a(`alphamid') d(`deltamid') n(`nunitmid') td(`nunitmid'-1);
powercal powernew, a(`alphamid') d(`deltamid') s(sdinfnew) n(`nunitmid') td(`nunitmid'-1);
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>0 if !missing(`Z');
};
assert powernew<=1 if !missing(powernew);
assert reldif(powernew,powerold)<=1e-9 if !missing(powernew);
* Note that sdinfnew has some extreme values *;
version 7:graph powerold sdinfnew, c(L) ylab yline(`powermid') xline(`sdinfmid') xlog
 t1("Calculations assuming Student's t-distribution");
more;
drop *new;

*
 Power and Nunit
*;
* Normal distribution, nunit to power (and back) *;
powercal powernew, a(`alphamid') d(`deltamid') s(`sdinfmid') n(nunitold);
powercal nunitnew, p(powernew) a(`alphamid') d(`deltamid') s(`sdinfmid') noceil;
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>=0 if !missing(`Z');
};
assert powernew<=1 if !missing(powernew);
assert reldif(nunitnew,nunitold)<1e-2 if !missing(nunitnew);
version 7:graph powernew nunitold, c(L) ylab xlab yline(`powermid') xline(`nunitmid')
 t1("Calculations assuming Normal distribution");
more;
drop *new;
* Normal distribution, power to nunit (and back) *;
powercal nunitnew, p(powerold) a(`alphamid') d(`deltamid') s(`sdinfmid') noceil;
powercal powernew, a(`alphamid') d(`deltamid') s(`sdinfmid') n(nunitnew);
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>=0 if !missing(`Z');
};
assert powernew<=1 if !missing(powernew);
assert reldif(powernew,powerold)<1e-9 if !missing(powernew);
version 7:graph powerold nunitnew, c(L) ylab xlab yline(`powermid') xline(`nunitmid')
 t1("Calculations assuming Normal distribution");
more;
drop *new;
* Student's t-distribution, nunit to power (and back) *;
powercal powernew, a(`alphamid') d(`deltamid') s(`sdinfmid') n(nunitold) td(nunitold-1);
powercal nunitnew, p(powernew) a(`alphamid') d(`deltamid') s(`sdinfmid') td(nunitold-1) noceil;
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>=0 if !missing(`Z');
};
assert powernew<=1 if !missing(powernew);
assert reldif(nunitnew,nunitold)<1e-2 if !missing(nunitnew);
version 7:graph powernew nunitold, c(L) ylab xlab yline(`powermid') xline(`nunitmid')
 t1("Calculations assuming Student's t-distribution");
more;
drop *new;

*
 Alpha and Delta
*;
* Normal distribution, delta to alpha (and back) *;
powercal alphanew, p(`powermid') d(deltaold) s(`sdinfmid') n(`nunitmid');
powercal deltanew, p(`powermid') a(alphanew) s(`sdinfmid') n(`nunitmid');
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>=0 if !missing(`Z');
};
assert alphanew<=1 if !missing(alphanew);
assert alphanew>=0 if !missing(alphanew);
assert reldif(deltanew,deltaold)<1e-9 if !missing(deltanew);
version 7:graph alphanew deltaold, c(L) ylab xlab yline(`alphamid') xline(`deltamid') ylog yreverse
 t1("Calculations assuming Normal distribution");
more;
drop *new;
* Normal distribution, alpha to delta (and back) *;
powercal deltanew, p(`powermid') a(alphaold) s(`sdinfmid') n(`nunitmid');
powercal alphanew, p(`powermid') d(deltanew) s(`sdinfmid') n(`nunitmid');
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>=0 if !missing(`Z');
};
assert alphanew<=1 if !missing(alphanew);
assert alphanew>=0 if !missing(alphanew);
assert reldif(alphanew,alphaold)<=1e-9 if !missing(alphanew);
version 7:graph alphaold deltanew, c(L) ylab xlab yline(`alphamid') xline(`deltamid') ylog yreverse
 t1("Calculations assuming Normal distribution");
more;
drop *new;
* Student's t-distribution, delta to alpha (and back) *;
powercal alphanew, p(`powermid') d(deltaold) s(`sdinfmid') n(`nunitmid') td(`nunitmid'-1);
powercal deltanew, p(`powermid') a(alphanew) s(`sdinfmid') n(`nunitmid') td(`nunitmid'-1);
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>=0 if !missing(`Z');
};
assert alphanew<=1 if !missing(alphanew);
assert alphanew>=0 if !missing(alphanew);
assert reldif(deltanew,deltaold)<1e-9 if !missing(deltanew);
version 7:graph alphanew deltaold, c(L) ylab xlab yline(`alphamid') xline(`deltamid') ylog yreverse
 t1("Calculations assuming Student's t-distribution");
more;
drop *new;
* Student's t-distribution, alpha to delta (and back) *;
powercal deltanew, p(`powermid') a(alphaold) s(`sdinfmid') n(`nunitmid') td(`nunitmid'-1);
powercal alphanew, p(`powermid') d(deltanew) s(`sdinfmid') n(`nunitmid') td(`nunitmid'-1);
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>=0 if !missing(`Z');
};
assert alphanew<=1 if !missing(alphanew);
assert alphanew>=0 if !missing(alphanew);
assert reldif(alphanew,alphaold)<=1e-9 if !missing(alphanew);
version 7:graph alphaold deltanew, c(L) ylab xlab yline(`alphamid') xline(`deltamid') ylog yreverse
 t1("Calculations assuming Student's t-distribution");
more;
drop *new;

*
 Alpha and Sdinf
*;
* Normal distribution, sdinf to alpha (and back) *;
powercal alphanew, p(`powermid') d(`deltamid') s(sdinfold) n(`nunitmid');
powercal sdinfnew, p(`powermid') a(alphanew) d(`deltamid') n(`nunitmid');
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>=0 if !missing(`Z');
};
assert alphanew<=1 if !missing(alphanew);
assert alphanew>=0 if !missing(alphanew);
assert reldif(sdinfnew,sdinfold)<1e-9 if !missing(sdinfnew);
version 7:graph alphanew sdinfold, c(L) ylab xlab yline(`alphamid') xline(`sdinfmid') ylog yreverse xlog
 t1("Calculations assuming Normal distribution");
more;
drop *new;
* Normal distribution, alpha to sdinf (and back) *;
powercal sdinfnew, p(`powermid') a(alphaold) d(`deltamid') n(`nunitmid');
powercal alphanew, p(`powermid') d(`deltamid') s(sdinfnew) n(`nunitmid');
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>=0 if !missing(`Z');
};
assert alphanew<=1 if !missing(alphanew);
assert alphanew>=0 if !missing(alphanew);
assert reldif(alphanew,alphaold)<=1e-9 if !missing(alphanew);
version 7:graph alphaold sdinfnew, c(L) ylab xlab yline(`alphamid') xline(`sdinfmid') ylog yreverse xlog
 t1("Calculations assuming Normal distribution");
more;
drop *new;
* Student's t-distribution, sdinf to alpha (and back) *;
powercal alphanew, p(`powermid') d(`deltamid') s(sdinfold) n(`nunitmid') td(`nunitmid'-1);
powercal sdinfnew, p(`powermid') a(alphanew) d(`deltamid') n(`nunitmid') td(`nunitmid'-1);
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>=0 if !missing(`Z');
};
assert alphanew<=1 if !missing(alphanew);
assert alphanew>=0 if !missing(alphanew);
assert reldif(sdinfnew,sdinfold)<1e-9 if !missing(sdinfnew);
version 7:graph alphanew sdinfold, c(L) ylab xlab yline(`alphamid') xline(`sdinfmid') ylog yreverse xlog
 t1("Calculations assuming Student's t-distribution");
more;
drop *new;
* Student's t-distribution, alpha to sdinf (and back) *;
powercal sdinfnew, p(`powermid') a(alphaold) d(`deltamid') n(`nunitmid') td(`nunitmid'-1);
powercal alphanew, p(`powermid') d(`deltamid') s(sdinfnew) n(`nunitmid') td(`nunitmid'-1);
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>=0 if !missing(`Z');
};
assert alphanew<=1 if !missing(alphanew);
assert alphanew>=0 if !missing(alphanew);
assert reldif(alphanew,alphaold)<=1e-9 if !missing(alphanew);
version 7:graph alphaold sdinfnew, c(L) ylab xlab yline(`alphamid') xline(`sdinfmid') ylog yreverse xlog
 t1("Calculations assuming Student's t-distribution");
more;
drop *new;

*
 Alpha and Nunit
*;
* Normal distribution, nunit to alpha (and back) *;
powercal alphanew, p(`powermid') d(`deltamid') s(`sdinfmid') n(nunitold);
powercal nunitnew, p(`powermid') a(alphanew) d(`deltamid') s(`sdinfmid') noceil;
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>=0 if !missing(`Z');
};
assert alphanew<=1 if !missing(alphanew);
assert alphanew>=0 if !missing(alphanew);
assert reldif(nunitnew,nunitold)<1e-9 if !missing(nunitnew);
version 7:graph alphanew nunitold, c(L) ylab xlab yline(`alphamid') xline(`nunitmid') ylog yreverse
 t1("Calculations assuming Normal distribution");
more;
drop *new;
* Normal distribution, alpha to nunit (and back) *;
powercal nunitnew, p(`powermid') a(alphaold) d(`deltamid') s(`sdinfmid') noceil;
powercal alphanew, p(`powermid') d(`deltamid') s(`sdinfmid') n(nunitnew);
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>=0 if !missing(`Z');
};
assert alphanew<=1 if !missing(alphanew);
assert alphanew>=0 if !missing(alphanew);
assert reldif(alphanew,alphaold)<=1e-9 if !missing(alphanew);
version 7:graph alphaold nunitnew, c(L) ylab xlab yline(`alphamid') xline(`nunitmid') ylog yreverse
 t1("Calculations assuming Normal distribution");
more;
drop *new;
* Student's t-distribution, nunit to alpha (and back) *;
powercal alphanew, p(`powermid') d(`deltamid') s(`sdinfmid') n(nunitold) td(nunitold-1);
powercal nunitnew, p(`powermid') a(alphanew) d(`deltamid') s(`sdinfmid') td(nunitold-1) noceil;
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>=0 if !missing(`Z');
};
assert alphanew<=1 if !missing(alphanew);
assert alphanew>=0 if !missing(alphanew);
assert reldif(nunitnew,nunitold)<1e-9 if !missing(nunitnew);
version 7:graph alphanew nunitold, c(L) ylab xlab yline(`alphamid') xline(`nunitmid') ylog yreverse
 t1("Calculations assuming Student's t-distribution");
more;
drop *new;

*
 Delta and Sdinf
*;
* Normal distribution, sdinf to delta (and back) *;
powercal deltanew, p(`powermid') a(`alphamid') s(sdinfold) n(`nunitmid');
powercal sdinfnew, p(`powermid') a(`alphamid') d(deltanew) n(`nunitmid');
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>=0 if !missing(`Z');
};
assert reldif(sdinfnew,sdinfold)<1e-9 if !missing(sdinfnew);
version 7:graph deltanew sdinfold, c(L) ylab xlab yline(`deltamid') xline(`sdinfmid') xlog
 t1("Calculations assuming Normal distribution");
more;
drop *new;
* Normal distribution, delta to sdinf (and back) *;
powercal sdinfnew, p(`powermid') a(`alphamid') d(deltaold) n(`nunitmid');
powercal deltanew, p(`powermid') a(`alphamid') s(sdinfnew) n(`nunitmid');
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>=0 if !missing(`Z');
};
assert reldif(deltanew,deltaold)<=1e-9 if !missing(deltanew);
version 7:graph deltaold sdinfnew, c(L) ylab xlab yline(`deltamid') xline(`sdinfmid') xlog
 t1("Calculations assuming Normal distribution");
more;
drop *new;
* Student's t-distribution, sdinf to delta (and back) *;
powercal deltanew, p(`powermid') a(`alphamid') s(sdinfold) n(`nunitmid') td(`nunitmid'-1);
powercal sdinfnew, p(`powermid') a(`alphamid') d(deltanew) n(`nunitmid') td(`nunitmid'-1);
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>=0 if !missing(`Z');
};
assert reldif(sdinfnew,sdinfold)<1e-9 if !missing(sdinfnew);
version 7:graph deltanew sdinfold, c(L) ylab xlab yline(`deltamid') xline(`sdinfmid') xlog
 t1("Calculations assuming Student's t-distribution");
more;
drop *new;
* Student's t-distribution, delta to sdinf (and back) *;
powercal sdinfnew, p(`powermid') a(`alphamid') d(deltaold) n(`nunitmid') td(`nunitmid'-1);
powercal deltanew, p(`powermid') a(`alphamid') s(sdinfnew) n(`nunitmid') td(`nunitmid'-1);
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>=0 if !missing(`Z');
};
assert reldif(deltanew,deltaold)<=1e-9 if !missing(deltanew);
version 7:graph deltaold sdinfnew, c(L) ylab xlab yline(`deltamid') xline(`sdinfmid') xlog
 t1("Calculations assuming Student's t-distribution");
more;
drop *new;

*
 Delta and Nunit
*;
* Normal distribution, nunit to delta (and back) *;
powercal deltanew, p(`powermid') a(`alphamid') s(`sdinfmid') n(nunitold);
powercal nunitnew, p(`powermid') a(`alphamid') d(deltanew) s(`sdinfmid') noceil;
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>=0 if !missing(`Z');
};
assert reldif(nunitnew,nunitold)<1e-9 if !missing(nunitnew);
version 7:graph deltanew nunitold, c(L) ylab xlab yline(`deltamid') xline(`nunitmid')
 t1("Calculations assuming Normal distribution");
more;
drop *new;
* Normal distribution, delta to nunit (and back) *;
powercal nunitnew, p(`powermid') a(`alphamid') d(deltaold) s(`sdinfmid') noceil;
powercal deltanew, p(`powermid') a(`alphamid') s(`sdinfmid') n(nunitnew);
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>=0 if !missing(`Z');
};
assert reldif(deltanew,deltaold)<=1e-9 if !missing(deltanew);
version 7:graph deltaold nunitnew, c(L) ylab xlab yline(`deltamid') xline(`nunitmid')
 t1("Calculations assuming Normal distribution");
more;
drop *new;
* Student's t-distribution, nunit to delta (and back) *;
powercal deltanew, p(`powermid') a(`alphamid') s(`sdinfmid') n(nunitold) td(nunitold-1);
powercal nunitnew, p(`powermid') a(`alphamid') d(deltanew) s(`sdinfmid') td(nunitold-1) noceil;
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>=0 if !missing(`Z');
};
assert reldif(nunitnew,nunitold)<1e-9 if !missing(nunitnew);
version 7:graph deltanew nunitold, c(L) ylab xlab yline(`deltamid') xline(`nunitmid')
 t1("Calculations assuming Student's t-distribution");
more;
drop *new;

*
 Sdinf and Nunit
*;
* Normal distribution, nunit to sdinf (and back) *;
powercal sdinfnew, p(`powermid') a(`alphamid') d(`deltamid') n(nunitold);
powercal nunitnew, p(`powermid') a(`alphamid') d(`deltamid') s(sdinfnew) noceil;
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>=0 if !missing(`Z');
};
assert reldif(nunitnew,nunitold)<1e-9 if !missing(nunitnew);
version 7:graph sdinfnew nunitold, c(L) ylab xlab yline(`sdinfmid') xline(`nunitmid') ylog
 t1("Calculations assuming Normal distribution");
more;
drop *new;
* Normal distribution, sdinf to nunit (and back) *;
powercal nunitnew, p(`powermid') a(`alphamid') d(`deltamid') s(sdinfold) noceil;
powercal sdinfnew, p(`powermid') a(`alphamid') d(`deltamid') n(nunitnew);
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>=0 if !missing(`Z');
};
assert reldif(sdinfnew,sdinfold)<=1e-9 if !missing(sdinfnew);
version 7:graph sdinfold nunitnew, c(L) ylab xlab yline(`sdinfmid') xline(`nunitmid') ylog
 t1("Calculations assuming Normal distribution");
more;
drop *new;
* Student's t-distribution, nunit to sdinf (and back) *;
powercal sdinfnew, p(`powermid') a(`alphamid') d(`deltamid') n(nunitold) td(nunitold-1);
powercal nunitnew, p(`powermid') a(`alphamid') d(`deltamid') s(sdinfnew) td(nunitold-1) noceil;
foreach Z of var *new {;
  summ `Z';assert !missing(r(N)) & (r(N)>0);
  assert `Z'>=0 if !missing(`Z');
};
assert reldif(nunitnew,nunitold)<1e-9 if !missing(nunitnew);
version 7:graph sdinfnew nunitold, c(L) ylab xlab yline(`sdinfmid') xline(`nunitmid') ylog
 t1("Calculations assuming Student's t-distribution");
more;
drop *new;


*
 Test sequence 5.
 Certify -powercal- against -sampsi-
 by calculating power and sample size both ways,
 calculating each of the parameters power and nunit
 as a function of each of the other 4 parameters.
*;

*
 Power from alpha
*;
gene double powernew=.;
local i1=0;
while `i1'<`nobsmax' {;
  local i1=`i1'+1;
  local alphacur=alphaold[`i1'];
  cap sampsi 0 `deltamid', onesample alpha(`alphacur') sd(`sdinfmid') n(`nunitmid');
  if _rc==0 {;
    qui replace powernew=r(power) in `i1';
  };
};
lab var powernew "power from sampsi";
powercal thingnew, a(alphaold) d(`deltamid') s(`sdinfmid') n(`nunitmid');
gene rdnew=reldif(powernew,thingnew);
summ rdnew;
version 7:graph thingnew powernew;
more;
assert rdnew<0.015 if !missing(rdnew);
drop *new;

*
 Power from delta
*;
gene double powernew=.;
local i1=0;
while `i1'<`nobsmax' {;
  local i1=`i1'+1;
  local deltacur=deltaold[`i1'];
  cap sampsi 0 `deltacur', onesample alpha(`alphamid') sd(`sdinfmid') n(`nunitmid');
  if _rc==0 {;
    qui replace powernew=r(power) in `i1';
  };
};
lab var powernew "power from sampsi";
powercal thingnew, a(`alphamid') d(deltaold) s(`sdinfmid') n(`nunitmid');
gene rdnew=reldif(powernew,thingnew);
summ rdnew;
version 7:graph thingnew powernew;
more;
assert rdnew<0.025 if !missing(rdnew);
drop *new;

*
 Power from sdinf
*;
gene double powernew=.;
local i1=0;
while `i1'<`nobsmax' {;
  local i1=`i1'+1;
  local sdinfcur=sdinfold[`i1'];
  cap sampsi 0 `deltamid', onesample alpha(`alphamid') sd(`sdinfcur') n(`nunitmid');
  if _rc==0 {;
    qui replace powernew=r(power) in `i1';
  };
};
lab var powernew "power from sampsi";
powercal thingnew, a(`alphamid') d(`deltamid') s(sdinfold) n(`nunitmid');
gene rdnew=reldif(powernew,thingnew);
summ rdnew;
version 7:graph thingnew powernew;
more;
assert rdnew<0.025 if !missing(rdnew);
drop *new;

*
 Power from nunit
*;
gene double powernew=.;
local i1=0;
while `i1'<`nobsmax' {;
  local i1=`i1'+1;
  local nunitcur=nunitold[`i1'];
  cap sampsi 0 `deltamid', onesample alpha(`alphamid') sd(`sdinfmid') n(`nunitcur');
  if _rc==0 {;
    qui replace powernew=r(power) in `i1';
  };
};
lab var powernew "power from sampsi";
powercal thingnew, a(`alphamid') d(`deltamid') s(`sdinfmid') n(nunitold);
gene rdnew=reldif(powernew,thingnew);
summ rdnew;
version 7:graph thingnew powernew;
more;
assert rdnew<0.015 if !missing(rdnew);
drop *new;

*
 Nunit from power
*;
gene double nunitnew=.;
local i1=0;
while `i1'<`nobsmax' {;
  local i1=`i1'+1;
  local powercur=powerold[`i1'];
  cap sampsi 0 `deltamid', onesample power(`powercur') alpha(`alphamid') sd(`sdinfmid');
  if _rc==0 {;
    qui replace nunitnew=r(N_1) in `i1';
  };
};
lab var nunitnew "nunit from sampsi";
powercal thingnew, p(powerold) a(`alphamid') d(`deltamid') s(`sdinfmid');
gene adnew=abs(nunitnew-thingnew);
summ adnew;
version 7:graph thingnew nunitnew;
more;
assert adnew<=1 if !missing(adnew);
drop *new;

*
 Nunit from alpha
*;
gene double nunitnew=.;
local i1=0;
while `i1'<`nobsmax' {;
  local i1=`i1'+1;
  local alphacur=alphaold[`i1'];
  cap sampsi 0 `deltamid', onesample power(`powermid') alpha(`alphacur') sd(`sdinfmid');
  if _rc==0 {;
    qui replace nunitnew=r(N_1) in `i1';
  };
};
lab var nunitnew "nunit from sampsi";
powercal thingnew, p(`powermid') a(alphaold) d(`deltamid') s(`sdinfmid');
gene adnew=abs(nunitnew-thingnew);
summ adnew;
version 7:graph thingnew nunitnew;
more;
assert adnew<=1 if !missing(adnew);
drop *new;

*
 Nunit from delta
*;
gene double nunitnew=.;
local i1=0;
while `i1'<`nobsmax' {;
  local i1=`i1'+1;
  local deltacur=deltaold[`i1'];
  cap sampsi 0 `deltacur', onesample power(`powermid') alpha(`alphamid') sd(`sdinfmid');
  if _rc==0 {;
    qui replace nunitnew=r(N_1) in `i1';
  };
};
lab var nunitnew "nunit from sampsi";
powercal thingnew, p(`powermid') a(`alphamid') d(deltaold) s(`sdinfmid');
gene adnew=abs(nunitnew-thingnew);
summ adnew;
version 7:graph thingnew nunitnew;
more;
assert adnew<=1 if !missing(adnew);
drop *new;

*
 Nunit from sdinf
*;
gene double nunitnew=.;
local i1=0;
while `i1'<`nobsmax' {;
  local i1=`i1'+1;
  local sdinfcur=sdinfold[`i1'];
  cap sampsi 0 `deltamid', onesample power(`powermid') alpha(`alphamid') sd(`sdinfcur');
  if _rc==0 {;
    qui replace nunitnew=r(N_1) in `i1';
  };
};
lab var nunitnew "nunit from sampsi";
powercal thingnew, p(`powermid') a(`alphamid') d(`deltamid') s(sdinfold);
gene adnew=abs(nunitnew-thingnew);
summ adnew;
version 7:graph thingnew nunitnew;
more;
assert adnew<=1 if !missing(adnew);
drop *new;


exit;