new; cls;
theta=2;       /* -1<=theta<=2 in this particular case */
lowerbound=0|0;
upperbound=1|1;
nn=14;                  /* max nn=14 in Light version of GAUSS */

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

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




/*
**  rej2method2.txt -  Bivariate random numbers generation from PDF given by acceptance-rejection method.
**                                                                             (Latin Hypercube version)
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**
**  Purpose:    Creates random numbers from any PDF given.
**
**  Format:     x=rej2method2(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.
*/


proc rej2method2(lowerbound,upperbound,nn,&f);
   local Points,s1,s2,fs,max,xx,i,y1,y2,U,UUU,f:proc;
   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);
/* 3D LHS data */
   UUU=LHss(nn,3);
/* acceptance-rejection method */
   xx=miss(0,0)~miss(0,0);
   i=1;
   do while i<=nn^3;
      y1=(upperbound[1]-lowerbound[1])*UUU[i,1]+lowerbound[1];
      y2=(upperbound[1]-lowerbound[1])*UUU[i,2]+lowerbound[1];
      U=UUU[i,3];
      if U<=f(y1,y2)/max;
         xx=xx|(y1~y2);
      endif;
      i=i+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;
   retp(xx);
endp;

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


/*
**  LHss.txt - Latin Hypercube U(0,1).
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Calculates Latin Hypercube numbers in a very easy way. new algorithm.
**
**  Format:     x=LHss(nn,dim);
**
**  Input:      nn       scalar,    max number of index  (1,2,3,...,nn)
**
**              dim      scalar,    dimension            (1,...,dim)
**
**
**  Output:     x        matrix ,   (nn^dim) x (dim) of resulting U(0,1) numbers
**
**  Notice:     It takes lots of memory to run. Light version of GAUSS will not work in most cases.
*/
proc LHss(nn,dim);
   local x,b,n,z,y,indexvec,d,U,i,index;
/* convert x(base 10) to y(base b) */
   x=seqa(1,1,nn^dim-1); 
   b=nn;
   n=maxc(log(x)/log(b))+1;
   z=reshape(b,n,rows(x));
   y=rev(recserrc(x,z))';
   indexvec=y+1;                     /* shift all elements by 1    */
   indexvec=ones(1,dim)|indexvec;    /* insert 1's at the 1st row  */
/* divide between 0 and 1 into nn segments and sample from each segment */
   d=1/nn;
   U=(d*indexvec-d*(indexvec-1)).*rndu(nn^dim,dim)+d*(indexvec-1);
   do while NOT(U>0);                                                      /* check if U contains 0's */
      U=(d*indexvec-d*(indexvec-1)).*rndu(nn^dim,dim)+d*(indexvec-1);
   endo;
/* randomize the order as a row */
   x=zeros(nn^dim,dim);
   i=1;
   do while i<=nn^dim-1;
      index=ceil((nn^dim-(i-1))*rndu(1,1));
      x[i,.]=U[index,.];
      if index==1;
         U=U[2:rows(U),.];
      elseif index==rows(U);
         U=U[1:rows(U)-1,.];
      else;
         U=U[1:(index-1) (index+1):rows(U),.];
      endif;
      i=i+1;
   endo;
   x[nn^dim,.]=U;
   retp(x);
endp;