proc iml;
options nodate pagesize=60 linesize=80;
*This program generates a sampling distribution of d values.
You need to input parameters mu (population means) SD (population SDs)
N (sample size) and number of samples Nsam;
*input mu, the population mean;
mu1 = 50;
mu2 = 55;
*input SD, the population standard deviation;
SD1= 10;
SD2= 20;
*Input the sample size per sample for the sampling distribution;
N1 = 25;
N2 = 50;
*Input the number of samples to be included in the empirical sampling distribution;
Nsam=1000;
* Derived parameters;
sigma = sqrt(((N1-1)#(SD1#SD1)+(N2-1)#(SD2#SD2))/(N1+N2-2));
delta = (mu1-mu2)/sigma;
print 'input' mu1 mu2 SD1 SD2 N1 N2 Nsam;
print 'pooled SD and population d' sigma delta;
*set up placeholder;
out1 = j(Nsam,1,-999);
*******************************************************************;
*  Here we go
*******************************************************************;
do h=1 to Nsam; * Set for Nsam iterations;
************************************;
*data generation for first sample;
************************************;
sample1 = normal(j(N1,1,0));
sample1 = sample1*SD1+j(N1,1,mu1);
mean1=sample1[+,]/N1;
V1=sample1-mean1;
V1 = j(N1,1,1)`*V1##2/(N1-1);
*****************************************;
* Data generation for second sample
*****************************************;
sample2 = normal(j(N2,1,0));
sample2 = sample2*SD2+j(N2,1,mu2);
mean2=sample2[+,]/N2;
V2=sample2-mean2;
V2 = j(N2,1,1)`*V2##2/(N2-1);
*****************************************;
*  Find pooled SD
*****************************************;
SDP = sqrt(((N1-1)#V1+(N2-1)#V2)/(N1+N2-2));
********************************************
*  Find d
********************************************;
meandiff = mean1-mean2;
d = meandiff/SDP;
*************************************
* Output
*************************************;
out1[h,1]=d;
end;
*create output datasets;
create out2 from out1;
append from out1;
quit;
data d2; set out2;
*If you want to save the distribution, write to an external file;
*file 'c:/sas/mike/corrs';
*put col1 8.4;
proc univariate plot normal;
run;