new; cls;
/* Custom Distribution (General version: End points need not to be zeros.)  */
/* Parameters to be changed */
xi={1,2,3,5};
yi={3,0,2,2};
nn=1000;

/* Calculate and Display the Result */
call customg(xi,yi,nn);
call customg2(xi,yi,nn);



/* composition method (faster) */

proc customg(xi,yi,nn);
   local n,a,b,i,h,na,x1,x2,x,suma,xx,index;
/* calculation of each area under the kinked lines. */
   n=rows(xi);
   a=zeros(n-1,2);
   b=zeros(n-1,1);
   i=1;
   do while i<=n-1;
      h=xi[i+1]-xi[i];
      if yi[i]<yi[i+1];
         a[i,1]=yi[i]*h;
         a[i,2]=(yi[i+1]-yi[i])*h/2;
         b[i]=1;
      elseif yi[i]>yi[i+1];
         a[i,1]=yi[i+1]*h;
         a[i,2]=(yi[i]-yi[i+1])*h/2;
         b[i]=2;
      else;
         a[i,1]=yi[i]*h;
         a[i,2]=0;
         b[i]=3;
      endif;
      i=i+1;
   endo;
   suma=sumc(sumc(a));
   a=a/suma; 
   na=round(a*nn);
   if sumc(sumc(na))/=nn;
      na[1,2]=na[1,2]+1;
   endif;
   x=miss(0,0);
   i=1;
   do while i<=n-1;
      if b[i]==1;
         if a[i,1]==0;
            x2=xi[i]+sqrt(rndu(na[i,2],1))*(xi[i+1]-xi[i]);
            x=x|x2;
         else;
            x1=(xi[i+1]-xi[i])*rndu(na[i,1],1)+xi[i];
            x2=xi[i]+sqrt(rndu(na[i,2],1))*(xi[i+1]-xi[i]);
            x=x|x1|x2;
         endif;
      elseif b[i]==2;
         if a[i,1]==0;
            x2=xi[i]+(1-sqrt(rndu(na[i,2],1)))*(xi[i+1]-xi[i]);
            x=x|x2;
         else;
            x1=(xi[i+1]-xi[i])*rndu(na[i,1],1)+xi[i];
            x2=xi[i]+(1-sqrt(rndu(na[i,2],1)))*(xi[i+1]-xi[i]); 
            x=x|x1|x2;
         endif;
      else;
         if a[i,1]/=0;
            x1=(xi[i+1]-xi[i])*rndu(na[i,1],1)+xi[i];
            x=x|x1;
         endif;
      endif;
      i=i+1;
   endo;
   x=x[2:rows(x)];
/* randomize the order (The same number could be drawn. Be careful.) */
   n=rows(x);
   xx=zeros(n,1);
   i=1;
   do while i<=n-1;
      index=ceil((n-(i-1))*rndu(1,1));
      xx[i]=x[index];
      if index==1;
         x=x[2:rows(x)];
      elseif index==rows(x);
         x=x[1:rows(x)-1];
      else;
         x=x[1:(index-1) (index+1):rows(x)];
      endif;
      i=i+1;
   endo;
   xx[n]=x;
   x=xx;
   x=x[1:nn];
/* descriptive statistics */
   print;
   print "   n=" rows(x);
   print "mean=" meanc(x);
   print "s.d.=" stdc(x);
   print " max=" maxc(x);
   print " min=" minc(x);
   print;
/* graph */
   library pgraph;
   graphset;
   pqgwin auto;
   title("Target Distribution");
      xy(xi,yi~(1/suma)*yi);
   title("Custom Relative Distribution");
      call hist(x,29);
   retp(x);
endp;



/* rejection method */

proc customg2(xi,yi,nn);
   local a,n,i,h,maxf,X,Y,U,min,max;
/* calculation of the area under the kinked lines. */
   a=0;
   n=rows(xi);
   i=1;
   do while i<=n-1;
      h=xi[i+1]-xi[i];
      a=a+(yi[i]+yi[i+1])*h/2;    /* sum of the areas of each trapezoid or triangle */
      i=i+1;
   endo;
/* rejection method */
   maxf=(1/a)*maxc(yi);
   min=xi[1];
   max=xi[n];
   X=zeros(nn,1);
   i=1;
   do while i<=nn;
      Y=(max-min)*rndu(1,1)+min;
      U=rndu(1,1);
      if U<=(1/a)*f(xi,yi,Y)/maxf;
         X[i]=Y;
         call sleep(0.1);
         print/lz i "out of " nn;
         i=i+1;
      endif;
   endo;
/* descriptive statistics */
   print;
   print "   n=" nn;
   print "mean=" meanc(X);
   print "s.d.=" stdc(X);
   print " max=" maxc(X);
   print " min=" minc(X);
   print;
/* graph */
   library pgraph;
   graphset;
   pqgwin auto;
   title("Target Distribution");
      xy(xi,yi~(1/a)*yi);
   title("Custom Relative Distribution");
      call hist(X,29);
   retp(X);
endp;



proc f(xi,yi,x);
   local index,y;
   index=sumc(xi.<x);
   y=(x-xi[index])/(xi[index+1]-xi[index])*(yi[index+1]-yi[index])+yi[index];
   retp(y);
endp;