new; cls;
a=2;
lowerbound=0;
upperbound=a;
n=10000;

fn f(x)=pdfsine(x,a);                 /* set PDF first (Sine Curve case) */
x=rejmethod(lowerbound,upperbound,n,&f);

library pgraph;
graphset;
pqgwin auto;
xy(seqa(0,a/1000,1001),pdfsine(seqa(0,a/1000,1001),a));
call hist(x,29);




/*
**  rejmethod.txt -  random numbers generation from PDF given by acceptance-rejection method.
** (C) Copyright 2002(Revised 2005) Yosuke Amijima. All Rights Reserved.
**
**
**  Purpose:    Creates random numbers from any PDF given.
**
**  Format:     x=rejmethod(lowerbound,upperbound,n,&f);
**
**  Input:      lowerbound  scalar, lower bound when p=0
**
**              upperbound  scalar, upper bound when p=1
**
**              n           scalar, number of random numbers
**
**              &f          a procedure, pdf function pre-setting as fn f(x)=pdf(x,para1,para2,...); 
**
**  Output:     x           n x 1 vector, random numbers
**
**
**  Notice:     Set 'f(x)=pdfname(x,para1,para2,....);' before using this procedure by regarding
**              parameter para1, para2... in original pdf function procedure as global variables.
*/

proc rejmethod(lowerbound,upperbound,n,&f);
   local Points,s,fs,max,x,i,y,U,f:proc;
   Points=10000;
/* get max f(x) roughly */
   s=seqa(lowerbound,(upperbound-lowerbound)/(Points-1),Points);
   fs=zeros(Points,1);
   i=1;
   do while i<=rows(fs);
      fs[i]=f(s[i]);
      i=i+1;
   endo;
   max=maxc(fs);
/* acceptance-rejection method */
   x=zeros(n,1);
   i=1;
   do while i<=n;
      y=(upperbound-lowerbound)*rndu(1,1)+lowerbound;
      U=rndu(1,1);
      if U<=f(y)/max; /* max=a*g(x) */
         x[i]=y;
         i=i+1;
      endif;
   endo;
   retp(x);
endp;

proc pdfsine(x,a);
   local d;
   d=pi/(2*a)*sin(pi*x/a);
   d=d.*((x.>=0).and (x.<=a));
   retp(d);
endp;