new; cls;
n=100;
mu={5 4};
Sig={1   0.9,
     0.9   1};
data=mu+rndn(n,2)*chol(Sig);
x=data[.,1]; y=data[.,2];
n=rows(x);
h0=1;      /* bandwidth for 2D kernel  */
points=99; /* # of points in 2D kernel */

{U1,U2}=Cnonpara(x,y,n,h0,points,"gauss",&F,&F);
library pgraph;
graphset;
_plctrl=-1;
xy(U1,U2);

proc F(x);
   local n,xi;
   n=rows(x);
   xi=rankindx(x,1);
   retp(xi/(n+1));
endp;


proc(2)=Cnonpara(x,y,nn,h0,points,types,&F1,&F2);
   local n,U,kern,point1,point2,k,i,j,U1,U2,F1:proc,F2:proc;
   n=rows(x);
   U=F1(x)~F2(y);
/* kernel Monte Carlo */
   kern=bikern(U,h0,points,types);
   kern=kern/maxc(maxc(kern));
   points=rows(kern);
   point1=seqa(0,1,points)/(points-1);
   point2=seqa(0,1,points)/(points-1);
   U=zeros(nn,2);
   k=1;
   do while k<=nn;
      i=rndu(1,1); j=rndu(1,1);
      if rndu(1,1)<kern[minindc(abs(point1-i)),minindc(abs(point2-j))];
         U[k,.]=i~j;
         k=k+1;
      endif;
   endo;
   U1=U[.,1]; U2=U[.,2];
   library pgraph;
   graphset;
   pqgwin auto;
   contour(point1',point2,kern);
   retp(U1,U2);
endp;

proc bikern(x,h0,points,types);
   local point1,point2,kern,i,j,x0; 
   point1=seqa(0,1,points)/(points-1);
   point2=seqa(0,1,points)/(points-1); 
   kern=zeros(points,points);
   i=1;
   do while i<=points;
      j=1;
      do while j<=points;
         x0=point1[i]~point2[j];
         kern[i,j]=kerneldens(x0,x,h0,types);
         j=j+1;
      endo;
      i=i+1;
   endo;
   retp(kern);
endp;

proc kerneldens(x0,x,h0,types);
   local n,hh,invs,dets,u,kv,kn;
   n=rows(x);
   hh=h0*n^(-1/6);
   invs=inv(chol(vcx(x)));
   dets=sqrt(det(vcx(x)));
   u=(x0-x)./hh;
   kv=kernels(u*invs,types);
   kn=prodc(kv');
   retp(sumc(kn)/(n*((hh^2)*dets)));
endp;

proc kernels(u,types);
   local kv;
   if types$=="uniform";
      kv=0.5*(abs(u).<=1); 
   elseif types$=="gauss";              
      kv=pdfn(u);
   elseif types$=="triangular";
      kv=(1-abs(u)).*(abs(u).<=1); 
   elseif types$=="biweight";         
      kv=(15/16)*((1-u^2)^2).*(abs(u).<=1);
   elseif types$=="triweight";
      kv=(35/32)*((1-u^2)^3).*(abs(u).<=1);
   elseif types$=="epanechnikov";
      kv=0.75*(1-u^2).*(abs(u).<=1);
   else;
      errorlog "ERROR:  Type must be either of the following:";
      errorlog "uniform/gauss/triangular/biweight/triweight/epanechnikov";
   endif;
   retp(kv);
endp;