           Kernel Density Estimation (Sheutrim Dec95)

/****Univariate and Bivariate Kernel Density Estimation Procedures***/
/********************************************************
                     Mainline code
          Sample Code to test the program below:
********************************************************/
library pgraph;

n = 10000;
x = 2*rndn(n,2);

let m     = 50 50;
let left  = -5 -5;
let right =  5  5;
let h     =  1  1;

bcnt = hist2(x,left,right,m);

fhat = dist2(bcnt,left,right,h,&kernele);

dg = (right-left)./(m-1);
g1 = seqa(left[1],dg[1],m[1]);
g2 = seqa(left[2],dg[2],m[2]);

surface(g2',g1,fhat);

bcnt = hist2(x[.,1],left[1],right[1],m[1]);

fhat = dist2(bcnt,left[1],right[1],h[1],&kerneln);

surface(g1,fhat);

end;

/************************** PROGRAM CODE *************************************/
/******************************************************************************
Univariate and Bivariate Kernel Density Estimation Procedures.
This code summarises the results of a request that I placed on Gaussians.  It
works well for my purposes and seems general enough to be of wider interest.
--Geoffrey Shuetrim

    kernelb
    kernelt
    kernelr
    kernele
    kerneln

    Dist1
    Hist1
    Dist2
    Hist2

These procedures have been compiled and modified by Geoffrey Shuetrim
on 12/3/95.  No responsibility for the code is taken by Geoffrey Shuetrim,
Berwin Turlach or Gary King.  This code is written and submitted for public
use.  There are no performance guarantees.Please acknowledge this code
(and its authors) if you find it useful in your own work.

Kernelb thru kerneln are modified versions of code developed by Gary King,
Department of Government, Harvard University. GKing@Harvard.Edu

Dist1, Dist2, Hist1 and Hist2 are modified versions of code developed by
B. Turlach. berwin@alphasun.anu.edu.au

Hist1 and Hist2 generate univariate and bivariate histograms using linear
binning methods.

 Inputs:     x     -- Nx1 or 2 matrix, the x values.
             left  -- the lower left bincenter.
             right -- the upper right bincenter.
             m     -- the number of bins in each direction.

 Output:     bincnts -- the histogram.

Dist1 and Dist2 generate univariate and bivariate density estimates using the
histograms generated by Hist1 and Hist2.  The particular kernel required must
be passed as a procedure name.

 Inputs:     bincnts -- the histogram.
             left    -- the lower left bincenter.
             right   -- the upper right bincenter.
             h       -- the window widths.
             @Kernel -- the particular kernel required.

 Output:     fhat    -- the density estimate.

*******************************************************************************/

@ BIWEIGHT KERNEL @
proc kernelb(z);
local a,res,t;
    t=(abs(z).<1);
    res=t.*((15/16)*((1-(z^2))^2));
    retp(res);
endp;

@ TRIANGULAR KERNEL @
proc kernelt(z);
local a,res,t;
    t=(abs(z).<1);
    res=t.*(1-abs(z));
    retp(res);
endp;

@ RECTANGULAR kernel @
proc kernelr(z);
local a,res,t;
    t=(abs(z).<1);
    res=t.*0.5;
    retp(res);
endp;

@ EPANECNIKOV KERNEL @
proc kernele(z);
local a,res,t;
    t=(abs(z).<sqrt(5));
    a=code(t,sqrt(5)|1);
    res=t.*((3/4)*(1-(1/5).*(z^2))./a);
    retp(res);
endp;

@ GAUSSIAN KERNEL @
proc kerneln(z);
local a,res,t;
    res = pdfn(z);
    retp(res);
endp;

PROC (1)=Dist1(bincnts,left,right,h,&Kernel);

/* Computes the matrix of kernel density estimate values from the
   binned bivariate data using window width h.*/

   LOCAL nbin,nobs,sh,kern,fhat,kernel:proc;

   nbin = ROWS(bincnts) ;
   nobs = SUMC(SUMC(bincnts)) ;
   sh = h * (nbin / (right - left)) ;

   kern = kernel(SEQA(1-nbin,1,2*nbin-1)/sh) ;

   fhat = CONV(bincnts,kern,nbin,2*nbin-1);

   fhat = fhat ./ (nobs * h);
   RETP(fhat) ;

ENDP;

PROC(1)=Hist1(x,left,right,M);
LOCAL step,
      range,
      ind,
      tind,
      bincnts,
      bcntr,
      w0,
      w1;

  range = right-left ; @ used as the distance to be analysed @

  @ Ensure that observations sitting at the extremes of the range are set
    equal to the boundaries of the range @
  x     = x - left' ;
  x     = x.*(x.>=0) - range' ;
  x     = x.*(x.<=0) + right' ;

  @ get the step size between bins @
  step  = range ./ (m-1) ;

  @ to take into account that we have two extra bins than grid points @
  m     = m+1 ;

  @ get bin indices for each observation
    ind is the index of the grid box that each observation lies in @
  ind   = floor((x-left')./step'+1) ;

  @ Get value of the upper right corner of the bin for each observation @
  bcntr = left' + ind .* step' ;

  @ Generate the weights to distribute frequency over bin centres @
  w00 = (bcntr - x) ./ step';

  @ Ensure that the weights are between zero and one @
  w0 = w0.*(w0.>=0) - 1;
  w0 = w0.*(w0.<=0) + 1;

@ Calculate the weights to be applied to each corner of the box @

  @ Weight on the lower left corner for each observation @
  w1 = 1-w0 ;

  @ Now construct the Histogram from the Data @

  @ Generate the breakpoints between the cell identifiers @
  tind = seqa(0,1,m);

  bincnts = countwts(ind,tind,w11);
  ind  = ind - 1;
  bincnts = bincnts + countwts(ind,tind,w0);

  RETP(bincnts[1:m-1]);

ENDP;

PROC (1)=Dist2(bincnts,left,right,h,&Kernel) ;

/* Computes the matrix of kernel density estimate values from the
   binned bivariate data using window widths h.*/

   LOCAL nbin1,nbin2,nobs,sh1,sh2,kern1,kern2,fhat,kernel:proc;

   nbin1 = ROWS(bincnts) ;
   nbin2 = COLS(bincnts) ;
   nobs = SUMC(SUMC(bincnts)) ;
   sh1 = h[1] * (nbin1 / (right[1,1] - left[1,1])) ;
   sh2 = h[2] * (nbin2 / (right[2,1] - left[2,1])) ;

   kern1 = kernel(SEQA(1-nbin1,1,2*nbin1-1)/sh1) ;
   kern2 = kernel(SEQA(1-nbin2,1,2*nbin2-1)/sh2) ;

   fhat = (CONV((CONV(bincnts',kern2,nbin2,2*nbin2-1))',
                               kern1,nbin1,2*nbin1-1))';

   fhat = fhat / (nobs * h[1] * h[2]) ;
   RETP(fhat) ;

ENDP;

PROC(1)=Hist2(x,left,right,M);
LOCAL step,
      range,
      ind,
      tind,
      bincnts,
      bcntr,
      w00,
      w01,
      w10,
      w11;

  range = right-left ; @ used as the distance to be analysed @

  @ Ensure that observations sitting at the extremes of the range are set
    equal to the boundaries of the range @
  x     = x - left' ;
  x     = x.*(x.>=0) - range' ;
  x     = x.*(x.<=0) + right' ;

  @ get the step size between bins @
  step  = range ./ (m-1) ;

  @ to take into account that we have two extra bins than grid points @
  m     = m+1 ;

  @ get bin indices for each observation
    ind is the index of the grid box that each observation lies in @
  ind   = floor((x-left')./step'+1) ;

  @ Get value of the upper right corner of the bin for each observation @
  bcntr = left' + ind .* step' ;

  @ Generate the weights to distribute frequency over bin centres @
  w00 = (bcntr - x) ./ step';

  @ Ensure that the weights are between zero and one @
  w00 = w00.*(w00.>=0) - 1;
  w00 = w00.*(w00.<=0) + 1;

@ Calculate the weights to be applied to each corner of the box @
  w10 = prodc(((1-w00[.,1])~w00[.,2])');
  w01 = prodc((w00[.,1]~(1-w00[.,2]))');
  w00 = prodc(w00') ;
  w11 = 1-w00-w01-w10 ;

  @ Now construct the Histogram from the Data @

  @ Specify the identifier for the bin containing each observation.
    We want to generate a vectorised histogram for the data so we use
    the fact that there is one row of cells for each column so multiply
    the row index by the number of columns and then add the column index.
    This number correctly identifies the cell containing the observation. @

  ind = ind[.,1]*m[2]+ind[.,2];

  @ Generate the breakpoints between the cell identifiers @
  tind = seqa(0,1,prodc(m));

  bincnts = countwts(ind,tind,w11);
  ind  = ind - 1;
  bincnts = bincnts + countwts(ind,tind,w10);
  ind  = ind - m[2];
  bincnts = bincnts + countwts(ind,tind,w00);
  ind  = ind + 1;
  bincnts = bincnts + countwts(ind,tind,w01);
  bincnts = reshape(bincnts,m[1],m[2]);

  RETP(bincnts[1:m[1]-1,1:m[2]-1]);

ENDP;
