/**********************************************************************/
/* Appendix A: Listing of Macro SERIALT */
/* Macro: SERIALT */
/* Description: Executes serial test on RANUNI */
/* Parameters: */
/* DIM # of dimensions (DIM-tuples) */
/* INT # of intervals along each dimension */
/* SAMPLES # of samples of random numbers */
/* NOFRV # of random vectors per sample */
/* SEED Seed of function RANUNI */
%macro serialt(dim=,int=,samples=,nofrv=,seed=0);
%do sample=1 %to &samples;
/* Generate the random vectors. File work has */
/* variables xi's (i=1,2,...,&dim), where xi */
/* is the interval number in which the random */
/* vector falls along the ith dimension. */
data work(drop=i);
do i=1 to &nofrv;
%do i=1 %to &dim;
x&i=1+int(&int*ranuni(&seed));
%end;
output;
end;
/* Count number of vectors in each hypercube */
proc freq data=work;
table x1 %if &dim>1 %then %do;
%do i=2 %to &dim;
* x&i %end; %end;
/ noprint out=work;
/* Calculate the serial test */
data work(keep=prob);
set work end=end;
retain chisq 0 expfreq;
if _n_=1 then expfreq=&nofrv/(&int**&dim);
chisq=chisq+(count-expfreq)**2/expfreq;
if end then do;
prob=1-probchi(chisq,&int**&dim-1);
output; end;
run;
/* Save the results of the serial test */
proc append base=results data=work;
%end;
proc format;
value signif 0.0000-0.0001='P<0.0001'
0.0001-0.001 ='P<0.001'
0.001 -0.01 ='P<0.01'
0.01 -0.05 ='P<0.05'
0.05 -0.10 ='P<0.10'
0.10 -0.15 ='P<0.15'
0.15 -1 ='P>0.15';
/* Summarize the results */
proc freq data=results;
table prob;
format prob signif.;
title1 "Result of the Serial Test";
title2 "-------------------------";
title3 "# of Samples: &samples";
title4 "# of Random Vectors/Sample: &nofrv";
title5 "# of Dimensions: &dim";
title6 "# of Intervals/Dimension: &int";
run;
%mend;
/**********************************************************************/
/* Appendix B: Listing of Macro TCOIN */
/* Macro: TCOIN */
/* Description: The macro simulates the */
/* coin toss and */
/* executes a special run test. */
/* Parameters: */
/* NTOSS Number of tosses */
/* P Probability of head */
/* SEED Seed of random number generator */
%macro tcoin(ntoss=,p=0.5,seed=0);
data outcome(keep=outcome);
array side(2) $1 sideh sidet;
retain sideh 'H' sidet 'T';
do i=1 to &ntoss;
outcome=side(rantbl(&seed,&p));
output;
end;
proc freq data=outcome;
table outcome;
/* Determine the runs */
data runs(keep=run);
set outcome;
by notsorted outcome;
retain run;
if first.outcome then run=0;
run=run+1;
if last.outcome then output;
/* Produce statistics of the runs */
proc univariate data=runs freq;
/* Produce run frequencies for the test */
proc freq data=runs;
table run / out=runfreq noprint;
/* Cut the run-distribution at 10 */
data runfreq;
set runfreq end=end;
retain rest 0;
if run<=10 then do; output;
return; end;
rest=rest+count;
if end then do; count=rest;
run=11;
output; end;
/* Calculate the expected run frequencies */
data expfreq;
retain cexpprob 0;
do run=1 to 10;
expprob=1/(2**run);
expfreq=(&ntoss/2)*expprob;
cexpprob=cexpprob+expprob;
output;
end;
run=11;
expfreq=(&ntoss/2)*(1-cexpprob);
output;
/* Calculate and print the goodness-of-fit test */
data _null_;
merge runfreq expfreq;
by run;
file print;
retain chisq 0;
if _n_=1 then
put 'RUN FREQ EXP.FREQ' / 24*'-';
if count=. then count=0;
chisq=chisq+(count-expfreq)**2/expfreq;
if run<=10
then put run 2. count 11. expfreq 11.;
else put '11+' count 10. expfreq 11.;
if run=11 then do;
prob=1-probchi(chisq,10);
put 24*'-' / 'RUN TEST:' chisq 15.4 /
'PROBABILITY:' prob 12.4; end;
%mend;
/**********************************************************************/
/* Appendix C: Listing of Macro TBCOIN */
/* Macro: TBCOIN */
/* Description: The macro simulates the */
/* "tossing of a biased coin" */
/* where the tosses are defined */
/* as two-coin flips */
/* Parameters: */
/* NTOSS Number of tosses */
/* P Probability of head */
/* SEED Seed of random number generator */
%macro tbcoin(ntoss=,p=,seed=0);
data outcome(keep=toss1 toss2);
array side(2) $1 sideh sidet;
retain sideh 'H' sidet 'T';
do i=1 to &ntoss;
toss1=side(rantbl(&seed,&p));
toss2=side(rantbl(&seed,&p));
if toss1^=toss2 then output;
end;
proc freq;
table toss1*toss2/list;
%mend;
/**********************************************************************/
/* Appendix D: Listing of Macro MHS1 */
/* Macro: MHS1 */
/* Description: Simulation of the Monty Hall */
/* problem */
/* Parameters: */
/* HOST Behavior of host. Possible values: */
/* SO - switch is always offered */
/* NSO - switch is never offered */
/* M - host is malevolent (switch is */
/* offered only when initial */
/* choice is a winning one) */
/* SPROB Probability of player to switch when */
/* switch is offered (0<=SPROB<=1) */
/* OUT output data set of the simulation */
/* NPLAYS Number of plays */
/* SEED Seed of function RANUNI */
%macro mhs1(host=so,sprob=1,out=out,nplays=,
seed=0);
data &out(drop=i);
array env(3) $1 env1-env3;
length kenv showenv initch finalch 3;
do i=1 to &nplays;
env(1)='E'; env(2)='E'; env(3)='E';
/* Place the key into an envelope */
kenv=rantbl(&seed,1/3,1/3);
env(kenv)='K';
/* Make initial choice */
initch=rantbl(&seed,1/3,1/3);
initprz =env(initch);
/* Host offers no switch */
showenv=.;
%if &host=NSO %then %do;
finalch=initch;
%end;
/* Host offers switch */
%if &host=SO %then %do;
link show;
link switch;
%end;
/* Host is malevolent */
%if &host=M %then %do;
if env(initch)='E'
then finalch=initch;
else do;
link show;
link switch;
end;
%end;
finalprz=env(finalch);
output;
end;
return;
%if &host=SO | &host=M %then %do;
/* Choose an empty envelope to show */
show:
if initch^=kenv
then showenv=6-(initch+kenv);
else do; showenv=rantbl(&seed,0.5);
if kenv=1 | (kenv=2 & showenv=2)
then showenv=showenv+1;
end;
return;
/* Make the switch */
switch:
if ranuni(&seed)<&sprob
then finalch=6-(initch+showenv);
else finalch=initch;
return;
%end;
run;
/* Present simulation results */
proc freq data=out;
table finalprz initprz*finalprz env1*env2*env3
env1*env2*env3*initch*initprz*showenv*
finalch*finalprz/list missing;
%mend;
/**********************************************************************/
/* Appendix E:
/* Macro MHS2 */
/* Description: Simulates a generalized Monty */
/* Hall problem */
/* Parameters */
/* NEMPTY number of empty envelopes */
/* NKEY number of envelops with key */
/* restriction: NEMPTY+NKEY<=200 */
/* NREMOVE number of empty envelopes to be removed */
/* NPLAYS number of games to simulate */
/* OFFER probability of host offering the */
/* switch */
/* SWITCH probability of player taking the */
/* switch */
/* SEED seed of the random number generator */
/* OUT name of the output dataset which */
/* has one variable (PROB) */
%macro mhs2(nempty=,nkey=,nremove=,nplays=,
offer=,switch=,seed=0,out=);
%let nenv=%eval(&nempty+&nkey);
data &out(keep=prob);
array env(&nenv) $1 env1-env&nenv;
length pack0 pack $&nenv;
/* Make the empty envelopes */
do i=1 to &nempty;
env(i)='E';
end;
/* Make the envelopes with key */
do i=&nempty+1 to &nenv;
env(i)='K';
end;
/* Shuffle the envelopes */
do i=1 to &nenv-1;
r=int((&nenv-i+1)*ranuni(&seed))+1;
temp=env(&nenv-i+1);
env(&nenv-i+1)=env(r);
env(r)=temp;
end;
/* pack0 saves the contents of all envelopes */
pack0=env(1);
do i=2 to &nenv;
pack0=trim(pack0)||env(i);
end;
/* Play the game &nplays times */
nwins=0;
do p=1 to &nplays;
pack=pack0;
/* Make initial selection */
sel=int(&nenv*ranuni(&seed))+1;
prize=substr(pack,sel,1);
if ranuni(&seed)<&offer and
ranuni(&seed)<&switch
then do;
/* Remove empty envelopes */
i=0; j=1;
do while (i<&nremove);
if j^=sel and env(j)='E'
then do;
i=i+1;
substr(pack,j,1)=' ';
end;
j=j+1;
end;
substr(pack,sel,1)=' ';
pack=compress(pack);
/* Make final selection */
sel=int((&nenv-&nremove-1)*ranuni(&seed))+1;
prize=substr(pack,sel,1);
end;
if prize='K' then nwins=nwins+1;
end;
/* Output the probability of winning */
prob=nwins/&nplays;
output;
run;
%mend;
/**********************************************************************/