new; cls;
theta=2;       /* -1<=theta<=2 in this particular case */
lowerbound=0|0;
upperbound=1|1;
n=1000;

fn f(u,v)=pdf(u,v,theta);
xx=rej2method(lowerbound,upperbound,n,&f);

library pgraph;
graphset;
pqgwin auto;
_plctrl=-1;
xy(xx[.,1],xx[.,2]);




/*
**  rej2method.txt -  Bivariate random numbers generation from PDF given by acceptance-rejection method.
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**
**  Purpose:    Creates random numbers from any PDF given.
**
**  Format:     x=rej2method(lowerbound,upperbound,n,&f);
**
**  Input:      lowerbound  2 x 1 vector, lower bounds when p=0
**
**              upperbound  2 x 1 vector, upper bounds 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 rej2method(lowerbound,upperbound,n,&f);
   local Points,s1,s2,fs,max,xx,i,y1,y2,U,f:proc;
   Points=1000;
/* get max f(x) roughly */
   s1=seqa(lowerbound[1],(upperbound[1]-lowerbound[1])/(Points-1),Points);
   s2=seqa(lowerbound[2],(upperbound[2]-lowerbound[2])/(Points-1),Points);
   fs=zeros(Points,1);
   i=1;
   do while i<=rows(fs);
      fs[i]=f(s1[i],s2[i]);
      i=i+1;
   endo;
   max=maxc(fs);
/* acceptance-rejection method */
   xx=zeros(n,2);
   i=1;
   do while i<=n;
      y1=(upperbound[1]-lowerbound[1])*rndu(1,1)+lowerbound[1];
      y2=(upperbound[2]-lowerbound[2])*rndu(1,1)+lowerbound[2];
      U=rndu(1,1);
      if U<=f(y1,y2)/max;
         xx[i,.]=y1~y2;
         i=i+1;
      endif;
   endo;
   retp(xx);
endp;

proc pdf(u,v,theta);
   retp(1+theta.*(6*u^2-6*u+1).*(6*v^2-6*v+1) );
endp;