data;
mux=50;
sigmax=10;
muy=50;
sigmay=10;
n=50;
rho=.50;
do i = 1 to n;
x=sigmax*rannor(0)+mux;
y=rho*((x-mux)/sigmax)+(sqrt(1-rho**2))*rannor(0); 
y=sigmay*y+muy;
output;
end;
proc gplot;
plot y*x;;
proc corr;
var x y;
run;
quit;