******************************************************* * MACRO: ONERAND * * ONE SAMPLE PITMAN RANDOMIZATION TEST FOR LOCATION. * * INPUT: dname, input data set name * * vfir, variable containing ovservations * * theta0, value of the median under H0 * * nrands, number of Monte Carlo samples * * OUTPUT: u, generalized Wilcoxon statistic * * pvu, pvd, pvb: upper, lower and 2-sided * * p-vals for generalized Wilcoxon test * * w, Wilcoxon statistic * * pvur, pvdr, pvbr: upper, lower and 2-sided * * p-vals for Wilcoxon test * *******************************************************; *options symbolgen mlogic mprint; %macro onerand(dname,vfir,theta0,nrands); * GET INITIAL VALUE OF TEST STATISTIC, W0; data datanew; set &dname; keep &vfir; run; data zdev; set datanew; devy=&vfir-&theta0; if devy<0 then devy=0; run; proc means data=zdev noprint; var devy; output out=zdevout sum=devsum; run; data _null_; set zdevout; call symput('w0',left(put(devsum,best10.))); run; data zabs; set datanew; absy=abs(&vfir-&theta0); sign=(sign(&vfir-&theta0)+1)/2; rabsy=absy; keep absy rabsy sign; run; proc rank data=zabs out=rabs; var rabsy; run; data rabs; set rabs; srabs=sign*rabsy; run; proc means data=rabs noprint; var srabs; output out=srabsout sum=rdevsum; run; data _null_; set srabsout; call symput('rw0',left(put(rdevsum,best10.))); run; proc means data=datanew noprint; var &vfir; output out=ndata n=n; run; data _null_; set ndata; call symput('nobs',left(put(n,best10.))); run; proc iml; *reset print; use zabs; read all var{absy} into z; use rabs; read all var{rabsy} into rz; pvu=0; pvd=0; pvb=0; pvur=0; pvdr=0; pvbr=0; do m=1 to &nrands; tempw=0; temprw=0; do j=1 to &nobs; l=uniform(0); if l>.5 then do; tempw=tempw+z[j,1]; temprw=temprw+rz[j,1]; end; end; if tempw<=&w0 then pvd=pvd+1; if tempw>=&w0 then pvu=pvu+1; if temprw<=&rw0 then pvdr=pvdr+1; if temprw>=&rw0 then pvur=pvur+1; end; pvu=pvu/&nrands; pvd=pvd/&nrands; pvur=pvur/&nrands; pvdr=pvdr/&nrands; x=pvd||pvu||pvdr||pvur; varnames={pvd,pvu,pvdr,pvur}; create zpvals from x [colname=varnames]; append from x; quit; data zzzout; set zpvals; u=&w0; w=&rw0; pvb=2*min(pvd,pvu); pvbr=2*min(pvdr,pvur); run; title 'Wilcoxon p-values'; proc print data=zzzout noobs; var w pvur pvdr pvbr; run; title 'Generalized Wilcoxon p-values'; proc print data=zzzout noobs; var u pvu pvd pvb; run; title ' '; %mend onerand; * This will do a million sample Monte Carlo; %onerand(ma2612.grind,diameter,.75,1000000);