new; cls;
theta=2;       /* -1<=theta<=2 in this particular case */
lowerbound=0|0;
upperbound=1|1;
nn=24;          /* Max nn depends on the shape of pdf. */

fn f(u,v)=pdf(u,v,theta);
xx=rej2method2s(lowerbound,upperbound,nn,&f);

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



/*
**  rej2method2s.txt -  Bivariate random numbers generation from PDF given by acceptance-rejection method.
**                                                               (Latin Hypercube version for GAUSS Light)
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**
**  Purpose:    Creates random numbers from any PDF given.
**
**  Format:     x=rej2method2s(lowerbound,upperbound,nn,&f);
**
**  Input:      lowerbound  2 x 1 vector, lower bounds when p=0
**
**              upperbound  2 x 1 vector, upper bounds when p=1
**
**              nn          scalar, number of partitions in each dimension
**
**              &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.
**
**              More bigger number of nn would be required in some inefficient case. The resulting
**              number of random number n may differ from time to time, and we cannot control it.
**
**              If you have GAUSS itself, you don't have to use this procedure. Only for GAUSS Light.
*/


proc rej2method2s(lowerbound,upperbound,nn,&f);
   local Points,s1,s2,fs,max,xx,i,y1,y2,U,UUU,f:proc;
   local d,j1,j2,j3,indexvec,nxx,index,x;
   Points=5000;
/* 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);
/* recursive LHS inside thre rejection method */
   xx=miss(0,0)~miss(0,0);
   d=1/nn;
   j1=0;
   do while j1<=nn-1;
      j2=0;
      do while j2<=nn-1;
         j3=0;
         do while j3<=nn-1;
            indexvec=(j1~j2~j3)+1;
               UUU=(d*indexvec-d*(indexvec-1)).*rndu(1,3)+d*(indexvec-1);
               y1=(upperbound[1]-lowerbound[1])*UUU[1]+lowerbound[1];
               y2=(upperbound[1]-lowerbound[1])*UUU[2]+lowerbound[1];
               U=UUU[3];
               if U<=f(y1,y2)/max;
                  xx=xx|(y1~y2);
               endif;
            j3=j3+1;
         endo;
         j2=j2+1;
      endo;
      j1=j1+1;
   endo;
   if rows(xx)==1;
      errorlog "ERROR: There is no point. Increase nn(The number of partition).";
      retp(".");
   else;
      xx=xx[2:rows(xx),.];
      print/lz "n=" rows(xx);
   endif;
/* randomize the order as a row */
   nxx=rows(xx);
   x=zeros(nxx,2);
   i=1;
   do while i<=nxx-1;
      index=ceil((nxx-(i-1))*rndu(1,1));
      x[i,.]=xx[index,.];
      if index==1;
         xx=xx[2:rows(xx),.];
      elseif index==rows(xx);
         xx=xx[1:rows(xx)-1,.];
      else;
         xx=xx[1:(index-1) (index+1):rows(xx),.];
      endif;
      i=i+1;
   endo;
   x[nxx,.]=xx;
   xx=x;
   retp(xx);
endp;

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