libname sassy "C:\CalPoly\Classes\418\materials\computing";

libname sassy "F:\_external\CalPoly\Classes\418\materials\computing";


************************************************************;
/* COVERAGE FOR WALD */

%macro cover_w(n,rep);

data sassy.cover_w_&n.;
n=&n;
rep=&rep;
do pi = 0.5, 0.45, 0.4, 0.35, 0.3, 0.25, 0.2, 0.15, 0.1, 0.05, 0.01;
    success=0;
    do i = 1 to rep;
        y = ranbin(0,n,pi);
        phat = y/n;
        low = 
        max(phat - 1.96*sqrt((phat*(1-phat))/n),0);
        upp = 
        min(phat + 1.96*sqrt((phat*(1-phat))/n),1);
        if low <= pi and pi <= upp then capture = 1;
        else capture = 0;
        success+capture;
    end;
    if pi = 0.5  then cov_&n._50_W = success/rep;
    if pi = 0.45 then cov_&n._45_W = success/rep;
    if pi = 0.4  then cov_&n._40_W = success/rep;
    if pi = 0.35 then cov_&n._35_W = success/rep;
    if pi = 0.3  then cov_&n._30_W = success/rep;
    if pi = 0.25 then cov_&n._25_W = success/rep;
    if pi = 0.2  then cov_&n._20_W = success/rep;
    if pi = 0.15 then cov_&n._15_W = success/rep;
    if pi = 0.1  then cov_&n._10_W = success/rep;
    if pi = 0.05 then cov_&n._05_W = success/rep;
    if pi = 0.01 then cov_&n._01_W = success/rep;
end;
run;

proc print data=sassy.cover_w_&n.;
run;

%mend;

%cover_w(100,1000000);
%cover_w(50,1000000);
%cover_w(25,1000000);
%cover_w(10,1000000);

************************************************************;
/* COVERAGE FOR SCORE */

%macro cover_s(n,rep);

data sassy.cover_s_&n.;
n=&n;
rep=&rep;
do pi = 0.5, 0.45, 0.4, 0.35, 0.3, 0.25, 0.2, 0.15, 0.1, 0.05, 0.01;
    success=0;
    do i = 1 to rep;
        y = ranbin(0,n,pi);
        phat = y/n;
        low = 
          max(
          (phat + 1.96**2/(2*n) - 
          1.96*sqrt(phat*(1-phat)/n + 1.96**2/(4*n**2)))/
          (1 + 1.96**2/n)
          ,0);

        upp = 
          min(
          (phat + 1.96**2/(2*n) + 
          1.96*sqrt(phat*(1-phat)/n + 1.96**2/(4*n**2)))/
          (1 + 1.96**2/n)
          ,1);
        if low <= pi and pi <= upp then capture = 1;
        else capture = 0;
        success+capture;
    end;
    if pi = 0.5  then cov_&n._50_S = success/rep;
    if pi = 0.45 then cov_&n._45_S = success/rep;
    if pi = 0.4  then cov_&n._40_S = success/rep;
    if pi = 0.35 then cov_&n._35_S = success/rep;
    if pi = 0.3  then cov_&n._30_S = success/rep;
    if pi = 0.25 then cov_&n._25_S = success/rep;
    if pi = 0.2  then cov_&n._20_S = success/rep;
    if pi = 0.15 then cov_&n._15_S = success/rep;
    if pi = 0.1  then cov_&n._10_S = success/rep;
    if pi = 0.05 then cov_&n._05_S = success/rep;
    if pi = 0.01 then cov_&n._01_S = success/rep;
end;
run;

proc print data=sassy.cover_s_&n.;
run;

%mend;

%cover_s(100,1000000);
%cover_s(50,1000000);
%cover_s(25,1000000);
%cover_s(10,1000000);

************************************************************;
/* COVERAGE FOR AGRESTI-COULL */

%macro cover_AC(n,rep);

data sassy.cover_AC_&n.;
n=&n;
rep=&rep;
do pi = 0.5, 0.45, 0.4, 0.35, 0.3, 0.25, 0.2, 0.15, 0.1, 0.05, 0.01;
    success=0;
    do i = 1 to rep;
        y = ranbin(0,n,pi);
        ptilde = (y+2)/(n+4);
        low = 
        max(ptilde - 1.96*sqrt((ptilde*(1-ptilde))/
           (n+4)),0);

        upp = 
        min(ptilde + 1.96*sqrt((ptilde*(1-ptilde))/
           (n+4)),1);

if low <= pi and pi <= upp then capture = 1;
        else capture = 0;
        success+capture;
    end;
    if pi = 0.5  then cov_&n._50_AC = success/rep;
    if pi = 0.45 then cov_&n._45_AC = success/rep;
    if pi = 0.4  then cov_&n._40_AC = success/rep;
    if pi = 0.35 then cov_&n._35_AC = success/rep;
    if pi = 0.3  then cov_&n._30_AC = success/rep;
    if pi = 0.25 then cov_&n._25_AC = success/rep;
    if pi = 0.2  then cov_&n._20_AC = success/rep;
    if pi = 0.15 then cov_&n._15_AC = success/rep;
    if pi = 0.1  then cov_&n._10_AC = success/rep;
    if pi = 0.05 then cov_&n._05_AC = success/rep;
    if pi = 0.01 then cov_&n._01_AC = success/rep;
end;
run;

proc print data=sassy.cover_AC_&n.;
run;

%mend;

%cover_AC(100,1000000);
%cover_AC(50,1000000);
%cover_AC(25,1000000);
%cover_AC(10,1000000);