/* Step 1 generate simulated data

SAS can't generate a Vector of beta varaites, though it can generate a vector of gamma.
So a vector of  beta(a,b) is either generated by Gamma(a,1)/(Gamma(a,1)+Gamma(b,1)) or using
rand("beta", a,b) in a do loop. Do loop is much faster in data step than IML*/

 
data a;
n=100000; a=3; b=2; 
c=PDF('beta',(a-1)/((a-1)+(b-1)),a,b); 
k=0;
do until (k=n);
  u=ranuni(0);
  x=ranuni(0);
    if (u<=PDF('beta',x,a,b)/c) then do;
      y=rand("beta",a,b); 
      k=k+1;
      output;;
    end;
end;
keep x y;
run;

/*step 2 kernel density estimate and plot of density*/

/* PROC KDE in SAS 9.1 or higher can generate estimate density and it plot at the same time
The one single Proc KDE will do the work as following:

ods html; ods graphics on;
proc kde data=new;
univar x y / plots = density;
run;
ods graphics off; ods html close; 

*/

/* for SAS 9.0 or lower. PROC KDE can't generate density plot, so we need Proc Gplot */

proc kde data=a out=one; var x; run;
proc kde data=a out=two; var y; run;
data new;
  merge one (rename=(density=x_density))
        two (rename=(density=y_density));
run;
goption hsize=5 vsize=5;
symbol c=green i=spline;
proc gplot data=new;
axis1   label=(h=1 'Simulated Beta') c=red;
axis2   label=(h=1 'True Beta') c=blue;
plot  x_density*x/haxis=axis1;
plot  y_density*y/haxis=axis2;
run;quit;