new; cls;
nu=1;
lowerbound=-4;
upperbound=4;
nn=70;         /* Try nn^2 points on 2D grid. */

fn f(x)=pdfged(x,nu);          /* set PDF first (Here, one-parameter GED case.) */
x=rejmethod2(lowerbound,upperbound,nn,&f);

library pgraph;
graphset;
pqgwin auto;
call hist(x,29);




/*
**  rejmethod2.txt -  random numbers generation from PDF given by acceptance-rejection method.
**                                                                   (Latin Hypercube version)
** (C) Copyright 2002(Revised 2005) Yosuke Amijima. All Rights Reserved.
**
**
**  Purpose:    Creates random numbers from any PDF given.
**
**  Format:     x=rejmethod2(lowerbound,upperbound,nn,&f);
**
**  Input:      lowerbound  scalar, lower bound when p=0
**
**              upperbound  scalar, upper bound 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 rejmethod2(lowerbound,upperbound,nn,&f);
   local Points,s,fs,max,x,i,y,U,UU,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);
/* 2D LHS data */
   UU=LH2ss(nn);
/* acceptance-rejection method */
   x=miss(0,0);
   i=1;
   do while i<=nn^2;
      y=(upperbound-lowerbound)*UU[i,1]+lowerbound;
      U=UU[i,2];
      if U<=f(y)/max;
         x=x|y;
      endif;
      i=i+1;
   endo;
   if rows(x)==1;
      errorlog "ERROR: There is no point. Increase nn(The number of partition).";
      retp(".");
   else;
      x=x[2:rows(x)];
      print/lz "n=" rows(x);
   endif;
   retp(x);
endp;

proc pdfged(x,nu);   /* This is pdf of one-parameter GED. */
   local lambda;
   if nu<=0;
      errorlog "ERROR: Parameter nu must be positive.";
      retp(".");
   endif;
   lambda=1/2^(1/nu).*(gamma(1/nu)./gamma(3/nu))^(1/2);
   retp( nu.*exp(-1/2*abs(x./lambda)^nu)./(lambda.*2^(1+1/nu).*gamma(1/nu)) );
endp;  



/*
**  LH2ss.txt - 2D 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=LH2ss(nn);
**
**  Input:      nn       scalar,    max number of index  (1,2,3,...,nn)
**
**  Output:     x        matrix ,   (nn^2) x 2 of resulting U(0,1) numbers
**
**  Notice:     It takes lots of memory to run.
*/
proc LH2ss(nn);
   local x,b,n,z,y,indexvec,d,U,i,index;
/* convert x(base 10) to y(base b) */
   x=seqa(1,1,nn^2-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,2)|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^2,2)+d*(indexvec-1);
   do while NOT(U>0);                                                      /* check if U contains 0's */
      U=(d*indexvec-d*(indexvec-1)).*rndu(nn^2,2)+d*(indexvec-1);
   endo;
/* randomize the order as a row */
   x=zeros(nn^2,2);
   i=1;
   do while i<=nn^2-1;
      index=ceil((nn^2-(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^2,.]=U;
   retp(x);
endp;