/* Simulating a joint probability */

/* Generate data */
data random;
call streaminit(12345);					/* Set the random number generator value */
do i = 1 to 1000;						/* Perform a loop from 1 to 1000 */
	y = rand("Normal",0,1);				/* Marginal distribution 1: std. normal */
	z = rand("Normal",5,2);				/* Marginal distribution 2: normal with mu=5 and sigma=2 */
	output;								/* Output to the data set "random" */
end;									/* Tell SAS that loop is finished */
drop i;									/* Drop the loop index, i */
run;

/* Calculate the joint probability of each outcome and plot the 
	bivariate density */
proc kde data=random;
bivar y z / out=bivar plots=(contour surface);
run;quit;


/* Check that the highest joint probability; it should occur at the
	means of each marginal distribution */
proc means data=bivar max;
var density;
run;