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}=Cempkern0(x,y,n,h0,points);
library pgraph;
graphset;
_plctrl=-1;
pqgwin auto;
xy(U1,U2);


proc(2)=Cempkern0(x,y,nn,h0,points);
   local n,xi,yj,index,U,kern,point1,point2,k,i,j,U1,U2;
   n=rows(x);
   xi=rankindx(x,1); yj=rankindx(y,1);
   index=xi~yj;
   U=index/n;
/* kernel Monte Carlo */
   kern=bikern(U,h0,points);
   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);
   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);
         j=j+1;
      endo;
      i=i+1;
   endo;
   retp(kern);
endp;

proc kerneldens(x0,x,h0);
   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=pdfn(u*invs); /* Gaussian kernel */
   kn=prodc(kv');
   retp(sumc(kn)/(n*((hh^2)*dets)));
endp;