/**************   UNIVARIATE KERNEL DENSITY ESTIMATORS  *************
In response to a request on gaussians for a kernel density
estimator, three members of the list volunteered code.
All three procedures a listed below, along with a few comments.
As usual, the code is for public, non-commercial use. No
performance guarantees are included. Please cite the authors
if you use the code.

The kernel density estimators posted by Gary King and David Baird
contain implementations of univarite kernel estimation via direct
implementation of the formulas defining the estimator.  Baird
comments that King's implementation is more complete than his
own, which was written for a specific project, but he nevertheless
graciously agree to archive it so that users would have access
to both implementations. Scott Thompson's procedure, which is
listed first, implements an alternative approach, using the fast
fourier transform. His experience is that this approach is
significantly faster than the direct implementation, especially when
the sample size is large and cross-validation of the density estimate
is performed to choose a smoothing parameter.
"The speedups for cross-validation are especially large.
The method exploits the fact that not much smoothing is needed when
there is a lot of data, in that it runs much faster when the bandwidth
parameter is very small than when it is large.

On the other hand, since the method involves a preliminary "histogram"
estimate, it introduces some approximation errors that could be large
when the output grid is coarse and/or the sample size is small."

All of the methods in this file are for univariate density
estimation. Thompson would like to be informed of any
generalizations of the fourier transform approach to
higher-dimensional density estimation problems, where
he expects it to be even more advantageous.
July 1996
*************************************************************/



/*************************************************************
UNIVARIATE KERNEL DENSITY ESTIMATOR:
 FAST FOURIER TRANSFORM APPROACH

Written by:
T. Scott Thompson            email:  thompson@atlas.socsci.umn.edu
Department of Economics              tsthomps@facstaff.wisc.edu
University of Minnesota      phone:  (608) 263-3867

The methods are described in Silverman's book, and in the Applied
Statistics articles cited in the source code below.  They exploit
the facts that a kernel smoother is a convolution of the kernel
with the empirical distribution of the data, that convolutions
are easy t manipulate via fourier transforms, that
cross-validation scores can be
approximated easily from the fourier transform of the density
estimate, and that the fourier transform of the Gaussian kernel has a
particularly simple form.

The FFT permits a partial work-around  for
facts that a kernel smoother is a convolution of the kernel with the
empirical distribution of the data, that convolutions are easy to
manipulate via fourier transforms, that cross-validation scores can be
the "curse of dimnesionality" that hits higher dimensional density
estimation.)

The code was inspired in part  by the comparable Fortran
routines in the Applied Statistics algorithm archive on Netlib.


Limitations of my code:

1) Univariate density estimation only.

2) The output must be computed on a fixed-width grid of points, where
the number of output points is a power of 2.  Finer grids produce
smaller approximation errors (unlike "direct" methods) but there do
not appear to be any real gains to computing output values at more
than about 512 points.  If output is desired at other points in the
same range then some kind of post-processing interpolation is probably
the way to go.

3) Only the Gaussian kernel is implemented.  The Gaussian kernel is
especially nice for this method since it has a closed form,
real-valued, smooth fourier transform function that dies off quickly
in the tails.  Generalization to other kernels with these properties
is trivial.  Generalization to symmetric kernels with
exponential-tailed FT is also easy, even when the FT must be computed
numerically.

Disclaimer: No warranty expressed or implied.  Use at your own risk!

Please cite me if you use this!  I put quite a bit of time into it.

********************************************************************/
/***************************************************************/
/*************  DENEST.SRC *************************************/
/***************************************************************/
/*
         Procedure __DEFFT              Go to line  31
         Procedure DENEST                          123
         Procedure DENESTCV                        251
*/

#define big 60    /* should be a number big enough so that exp(-big)
                     is effectively zero but small enough so that
                     calculating exp(-big) will not cause underflow */

declare matrix __dataft ?= 0;






/*---------------------------------------------------------------------*/
/*      PROC __DEFFT   Calculate FFT of discretized univariate data   */
/*---------------------------------------------------------------------*/
/*

Note:  This is a utility routine to do the initial setup for DENEST
       and DENESTCV.  The discretization scheme is taken from:

             M.C. Jones and H.W. Lotwick, "Remark AS R50.  A Remark on
             Algorithm AS 176.  Kernel Density Estimation Using the Fast
             Fourier Transform," Applied Statistics 33 (1984), 120-122.

Usage:  __defft( dt, dlo, dhi, nft );

Input:  DT        = N x 1  vector of input data points.

        DLO, DHI  = Scalars:  low and high values of interval on which the
                             smooth is to be computed.

        WINDOW    = Smoothing parameter

        NFT       = Number of output points desired.  Must be a power
                    of two between 32 and 65,536.

Output:  The global __dataft is (re)computed.

Globals: __dataft = NFT x 2 vector containing the Fourier transform of
                    the data (after discretizing).  The purpose of this
                    routine is to initialize this global.

Important:  __DEFFT must be called whenever DATA, DLO, DHI, or NFT
            are changed since subsequent calls to the other routines
            will produce invalid results otherwise.

*/

proc (0) = __defft( dt, dlo, dhi, nft );

   local ndt, step, i, dlo1, bin, k, err13, ftr, fti;

   if (nft < 32) or (nft > 65536);
      errorlog("Error in __DEFFT: nft must be between 32 and 65,536.");
      end;
   endif;

   nft = maxc( nft | .5 );
   i = nft;
   do while i > 1;
      i = i / 2;
   endo;
   if i < 1;
      errorlog("Error in __DEFFT: NFT must be a power of 2.");
      end;
   endif;

   if not (dhi > dlo);
      errorlog("Error in _DEFFT: DHI should be larger than DLO.");
      end;
   endif;

   step = (dhi - dlo) / nft;           /*  step size for output */
   ndt = rows(dt);
   __dataft = zeros(nft+1,1);
   dlo1 = dlo - step/2;
   dt = (dt - dlo1) / step;
   bin = trunc(dt);
   dt = dt - bin;
   bin = bin | (bin + 1);
   dt = ( (1-dt) | dt ) / (ndt*step);
   err13 = error(13);
   i = 1;
   do until i == nft+1;
      k = indexcat(bin,i);
      if not scalmiss(k);
         __dataft[i] = sumc( dt[k] );
      endif;
      i = i+1;
   endo;
   k = indexcat(bin,nft);
   if not scalmiss(k);
      __dataft[1] = __dataft[1] + sumc( dt[k] );
   endif;
   __dataft = __dataft[1:nft];
   clear dt, bin, k;
   {ftr,fti} = rfft(__dataft);  /* compute discrete Fourier transform
                                   of the data */
   __dataft = ftr ~ fti;
   endp;




/*------------------------------------------------------------------------*/
/*      PROC DENEST     Estimate density using kernel smoother via FFT    */
/*------------------------------------------------------------------------*/
/*
Usage:  smooth = denest( data, dlo, dhi, window, nft );

Input:  DATA      = N x 1  vector of input data points.

        DLO, DHI  = Scalars:  low and high values of interval on which the
                             smooth is to be computed.

        WINDOW    = scalar.  Smoothing parameters

        NFT       = Number of output points desired.  Must be a power
                    of two between 32 and 65,536.

Output: SMOOTH    = NFT x 1 vector of density estimates with rows
                    corresponding to the evenly spaced points given by:

                         DLO + (I-1/2)*(DHI-DLO)/NFT, I=1,...,NFT

Globals: __dataft = NFT x 2 vector containing the Fourier transform of
                    the data (after discretizing).

Important:  __dataft must be re-initialized to a scalar zero whenever
            DATA, DLO, DHI, or NFT are changed since subsequent calls will
            produce invalid results otherwise.  Alternatively, the user may
            call __defft directly to reinitialize the array.

Note:  An alternative method for forcing initialization of __dataft
       is to replace NFT with the negative of the desired value.

Adapted from Applied Statistics Algorithm AS 176 as modified by R 50.

References:  B.W. Silverman, "Algorithm 176.  Kernel Density Estimation using
             the Fast Fourier Transform," Applied Statistics 31 (1982), 93-99.

             M.C. Jones and H.W. Lotwick, "Remark AS R50.  A Remark on
             Algorithm AS 176.  Kernel Density Estimation Using the Fast
             Fourier Transform," Applied Statistics 33 (1984), 120-122.

*/

proc (1) = denest( data, dlo, dhi, window, nft );

   local i, step, smooth, nft2, hw, fac, fac1, jhi, jmax;

   if window LE 0;
      errorlog("Error in DENEST: Window width must be positive.");
      end;
   endif;

   if nft < 0;                  /* force reinitialization if nft < 0 */
      __dataft = zeros(1,1);
      nft = abs(nft);
   endif;

   if (nft < 32) or (nft > 65536);
      errorlog("Error in DENEST: nft must be between 32 and 65,536.");
      end;
   endif;

   nft = maxc( nft | .5 );
   i = nft;
   do while i > 1;
      i = i / 2;
   endo;
   if i < 1;
      errorlog("Error in DENEST: NFT must be a power of 2.");
      end;
   endif;

   if not (dhi > dlo);
      errorlog("Error in DENEST: DHI should be larger than DLO.");
      end;
   endif;

/*  If __dataft is scalar then discretize data and initialize __dataft.  */

   if rows(__dataft)*cols(__dataft) == 1;

      __defft( data, dlo, dhi, nft );

   else;               /* check to see if __dataft has proper dimensions */

      if not (nft == rows(__dataft));
         errorlog("Error in DENEST: NFT should match rows(__dataft).");
         end;
      endif;

      if cols(__dataft) /= 2;
         errorlog("Error in DENEST: __dataft should have two columns.");
         end;
      endif;

   endif;

   step = (dhi - dlo) / nft;         /*  step size for output */
   nft2 = nft / 2;                   /*  nft2 is an integer           */
   hw = window / step;
   fac1 = pi * hw / nft;             /*   hw/nft = window/(dhi-dlo)   */
   fac1 = 2 * fac1 * fac1;           /*   (1/2) * (2*pi*hw/nft)**2    */

/*     Now find FFT of the density estimates */

   smooth = zeros(nft,2);
   jhi = trunc(sqrt( big / fac1 ));
   jmax = minc( nft2-1 | jhi );
   smooth[1,.] = __dataft[1,.];
   fac = seqa(1,1,jmax);
   fac = fac .* fac * fac1;
   fac = exp( -fac) ./ ( 1 - fac / (6*hw*hw) );
   smooth[2:jmax+1,.] = __dataft[2:jmax+1,.] .* fac;
   if jhi > nft2-1;
      smooth[nft2+1,1] = exp( -fac1 * nft2*nft2 ) * __dataft[nft2+1,1];
   endif;
   smooth[nft2+2:nft,.] = rev(smooth[2:nft2,1]~-smooth[2:nft2,2]);

/*    Invert Fourier transform of SMOOTH to get estimate  */

   smooth = rffti(smooth[.,1],smooth[.,2]);
   retp( smooth .* (smooth .GE 0) );          /* zero negative values. */

   endp;




/*------------------------------------------------------------------------*/
/*      PROC DENESTCV   Calculate kernel density CV score.                */
/*------------------------------------------------------------------------*/
/*
Usage:  score = denestcv( data, dlo, dhi, window, nft );

Input:  DATA      = N x 1  vector of input data points.

        DLO, DHI  = Scalars:  low and high values of interval on which the
                             smooth is to be computed.

        WINDOW    = J x 1 vector of smoothing parameters

        NFT       = Number of output points desired.  Must be a power
                    of two between 32 and 65,536.

Output: SCORE     = J x 1 vector of cross-validation scores corresponding
                    to the various elements of WINDOW.  The scores are
                    calculated according to a corrected form of formula
                    (3.62), p. 65 of Silverman (1986).

                    The score will be minimized at values of WINDOW that
                    minimize an approximation to

                                Integral (fhat - f)^2

                    where fhat is the kernel density function (as would be
                    estimated by DENEST) and f is the true density.

Globals: __dataft = NFT x 2 vector containing the Fourier transform of
                    the data (after discretizing).

Important:  See the discussion of __dataft in procedure DENEST.

Notes:  1)  This routine does not incorporation the bias correction
            factor incorporated into routine DENEST.

        2)  There is an error in Silverman's equations 3.61 and 3.62
            which has been corrected here.

References:  B.W. Silverman (1986).  Density Estimation for Statistics
             and Data Analysis.  (Chapman and Hall: London), pp. 61-66.

*/

proc (1) = denestcv( data, dlo, dhi, window, nft );

   local i, step, ndt, nw, nft2, hw, fac, fac1, jhi1, jhi2, jmax, cv;

   if minc(window) LE 0;
      errorlog("Error in DENESTCV: Window width must be positive.");
      end;
   endif;

   if nft < 0;                  /* force reinitialization if nft < 0 */
      __dataft = zeros(1,1);
      nft = abs(nft);
   endif;

   if (nft < 32) or (nft > 65536);
      errorlog("Error in DENESTCV: nft must be between 32 and 65,536.");
      end;
   endif;

   nft = maxc( nft | .5 );
   i = nft;
   do while i > 1;
      i = i / 2;
   endo;
   if i < 1;
      errorlog("Error in DENESTCV: NFT must be a power of 2.");
      end;
   endif;

   if not (dhi > dlo);
      errorlog("Error in DENESTCV: DHI should be larger than DLO.");
      end;
   endif;

/*  If __dataft is scalar then discretize data and initialize __dataft.  */

   if rows(__dataft)*cols(__dataft) == 1;

      __defft( data, dlo, dhi, nft );

   else;               /* check to see if __dataft has proper dimensions */

      if not (nft == rows(__dataft));
         errorlog("Error in DENESTCV: NFT should match rows(__dataft).");
         end;
      endif;

      if cols(__dataft) /= 2;
         errorlog("Error in DENESTCV: __dataft should have two columns.");
         end;
      endif;

   endif;

   ndt = rows(data);
   nw = rows(window);
   step = (dhi - dlo) / nft;         /*  step size for output */
   nft2 = nft / 2;                   /*  nft2 is an integer           */
   hw = window ./ step;
   fac1 = pi * hw ./ nft;            /*   hw/nft = window/(dhi-dlo)   */
   fac1 = 2 * fac1 .* fac1;          /*   (1/2) * (2*pi*hw/nft)**2    */

/*  Now find the CV score */

   jhi1 = minc( trunc(sqrt( big ./ fac1'     )) |
                nft2.*ones(1,nw)                );
   jhi2 = minc( trunc(sqrt( big ./ (2*fac1') )) |
                nft2.*ones(1,nw)                );
   jmax = maxc(maxc( jhi1' | jhi2' ));
   cv = zeros(jmax,nw);
   fac = seqa(1,1,jmax);
   fac = fac .* fac .* fac1';
   i = 1;
   do while i LE nw;
      cv[1:jhi2[i],i] =                     exp(-2*fac[1:jhi2[i],i]);
      cv[1:jhi1[i],i] = cv[1:jhi1[i],i] - 2*exp(-  fac[1:jhi1[i],i]);
      i = i+1;
   endo;
   cv = cv .* ( (__dataft[2:jmax+1,.].*__dataft[2:jmax+1,.])*(1|1) );
   retp( 2.*(dhi-dlo).*sumc(cv) - 1./(dhi-dlo)
                                + 2./(ndt.*window.*sqrt(2.*pi)) );
   endp;




/**************************************************************
** DENS.G - KERNEL DENSITY ESTIMATION
** Gary King, Department of Government, Harvard University
** <GKing@Harvard.Edu>
**
**   Format:    {px,py}=dens(y);
**
**   input:     y = Nx1 vector
**
**   output:    plot of y (horizontally) by its est'd density f(y) (vertically).
**              Whiskers are also plotted along horizontal axis for each datum
**
**              px = x coordinates for the plot (_pts x 1 vector)
**              py = y coordinates for the plot (_pts x 1 vector)
**
**    Global values are as follows, and may be changed:
**   _smth=0;         smoothing parameter>0. the larger, the more smoothing.
**                    =0 for automatic calculation of density.
**                    =-1 for automatic calculation for description.
**
**   _strt=0;         x axis minimum. (if _strt=_endd, set at min and max).
**   _endd=0;         x axis maximum.
**   _pts=100;        number of points to plot
**   _whiskr=-1;      -1 = draw whiskers at 1/10th size;
**                    >0 = draw whiskers at size _whiskr;
**                     0 = do not draw whiskers
**                    Whiskers are always drawn at y==0, so you may have to
**                    use ytics to scale the y-axis.
**   _jitter=0;       0=nothing extra. #=add # jitter to each
**                      whisker if _whiskr/=0;.
**
**   let _kern=E;     kernel type. N=normal, E=Epanechnikov,
**                                 B=biweight, T=triangular, R=rectangular
**                                 TN=Doubly Truncated Normal
**  for option _kern=="TN",
**      _Tleft=_strt;   truncate distribution on left at _Tleft
**      _Tright=_endd;  truncate distribution on right at _Tright
**
**   _output = 1;     1 = print density plot (default); 0 = do not print.
**
**  Reference:  B.W. Silverman. 1986. _Density Estimation for Statistics
**      and Data Analysis_.  London: Chapman and Hall.
**
**  An Application: Gary King. "Constituency Service and Incumbency
**      Advantage,"  _British Journal of Political Science_, 21, 1
**      (January, 1991): 119--128.
********************************************************************/
declare matrix _smth ?= 0;
declare matrix _strt ?= 0;
declare matrix _endd ?= 0;
declare matrix _pts  ?= 100;
declare matrix _Tleft ?= 0;
declare matrix _Tright ?= 0;
declare string _kern ?= "E";
declare matrix _whiskr ?= -1;
declare matrix _jitter ?= 0;
declare matrix _output ?= 1;
external proc code,xy,seqas;
external matrix _pline;
@ Kernels @
@ arguments of kernels: z=input value, m=mean, h=smoothing parameter _smth @
proc kerneln(z,m,h);local res;z=(z-m)./h;      @ NORMAL kernel @
    res=(1./sqrt(2*pi))*exp(-(1/2).*(z^2));retp(res);endp;
proc kernelTN(z,m,h,Tleft,Tright);
    local t,tl,res,tr,zz;zz=(z-m)./h;          @ TRUNCATED NORMAL kernel @
    tl=(Tleft-m)./h; tr=(Tright-m)./h;
    t=((z.>Tleft).and (z.<Tright));
    res=t.*pdfn(zz)./(1-cdfn(tl)-(1-cdfn(tr)));
    retp(res);endp;
proc kernele(z,m,h);local a,res,t;z=(z-m)./h;  @ EPANECHNIKOV kernel @
    t=(abs(z).<sqrt(5)); a=code(t,sqrt(5)|1);
    res=t.*((3/4)*(1-(1/5).*(z^2))./a);retp(res);endp;
proc kernelb(z,m,h);local a,res,t;z=(z-m)./h;  @ BIWEIGHT kernel @
    t=(abs(z).<1);
    res=t.*((15/16)*((1-(z^2))^2));retp(res);endp;
proc kernelt(z,m,h);local a,res,t;z=(z-m)./h;  @ TRIANGULAR kernel @
    t=(abs(z).<1);
    res=t.*(1-abs(z));retp(res);endp;
proc kernelr(z,m,h);local a,res,t;z=(z-m)./h;  @ RECTANGULAR kernel @
    t=(abs(z).<1);
    res=t.*0.5;retp(res);endp;

@ creates two nx1 vectors to plot @
proc (2)=density(y,strt,endd,pts,h,kerna,Tleft,Tright);
    local i,t,px,py;    local kern:proc;
    kerna=upper(kerna);
    px=seqas(strt,endd,pts);    py=px;
    format /rdn 4,0;
    "Kernel Density Estimation";
    "Calculating ";;pts;;" Points: ";;
    i=1;do while i<=pts;
        if     (kerna$=="N"); t=kernelN(py[i,1],y,h)./h;
        elseif (kerna$=="E"); t=kernelE(py[i,1],y,h)./h;
        elseif (kerna$=="B"); t=kernelB(py[i,1],y,h)./h;
        elseif (kerna$=="T"); t=kernelT(py[i,1],y,h)./h;
        elseif (kerna$=="R"); t=kernelR(py[i,1],y,h)./h;
        elseif (kerna$=="TN");t=kernelTN(py[i,1],y,h,Tleft,Tright)./h;
        else; errorlog "_kernel specified incorrectly";stop;endif;
        py[i,1]=sumc( t )/rows(y);
        if i==10*floor(i/10);i;;endif;
    i=i+1;endo;?;retp(px,py);endp;

@ Kernel Density Estimate; sets defaults, calls density, then xy @
proc 2=dens(y);local strt,endd,h,kern,pts,px,py,smth,std,Tleft,Tright,os;
    if cols(y)/=1; errorlog "Argument must be a column vector";end;endif;
    strt=_strt;endd=_endd;
    if     strt>endd;errorlog "error: _strt>_endd";end;
    elseif strt==endd;strt=minc(y);endd=maxc(y);endif;
    pts=int(_pts);
    if pts<=2;errorlog "_pts must be greater than 2.  Try _pts=100;";endif;
    smth=_smth;
    if (smth<0).and(smth/=-1);errorlog "_smth must be -1 or > than 0";stop;
    elseif smth==0; py=sortc(y,1);
           std=minc((py[int(3*rows(py)/4)]-py[int(rows(py)/4)])/1.34|stdc(py));
           smth=0.9*std*(rows(py)^(-0.2));
    elseif smth==-1; py=sortc(y,1);
           std=minc((py[int(3*rows(py)/4)]-py[int(rows(py)/4)])/1.34|stdc(py));
           smth=0.9*std*(rows(py)^(-0.2))/2; endif;
    kern=_kern;
    Tleft=0;Tright=0;
    if kern$=="TN";
        if _Tleft==_Tright;  Tleft=strt;   Tright=endd;
        else;                Tleft=_Tleft; Tright=_Tright; endif;
    endif;
    {px,py}=density(y,strt,endd,pts,smth,kern,Tleft,Tright);
    std=rndu(rows(y),1)*_jitter-(_jitter/2);
    y=y+std;
    if _output==1;
    if _whiskr==-1;
        os=ones(rows(y),1);
        _pline=os~ (os*6)~ y~ (os*0)~ y~ (os*(maxc(py)/15))~ os~ (os*15)~(os*0);
        clear os;
    elseif _whiskr>0;
        os=ones(rows(y),1);
        _pline=os~ (os*6)~ y~ (os*0)~ y~ (os*_whiskr)~ os~ (os*15)~(os*0);
        clear os;
    endif;
        xy(px,py);
        format/rd 5,4;
        if rows(_smth)==1;?;"Smoothing Parameter: _smth=";;smth;endif;
    endif;
    retp(px,py);endp;

/*
**  create vector of PTS evenly spaced points between STRT and ENDD,
**  not including the end points.
*/
proc seqas(strt,endd,pts);local t;t=(endd-strt)/pts;
    retp( seqa(strt+0.5*t,t,pts) );endp;






/************************************************************
**************************************************************
UNIVARIATE KERNEL DENSITY ESTIMATOR

  Dr David Baird  Biometrician
  Mail:           AgResearch, PO Box 60, Lincoln, NEW ZEALAND
  Phone:          +64 3 325 6900    Fax: +64 3 325 2946
  Direct Dial:    +64 3 325 6901 3975#
  Internet:       BairdD@AgResearch.CRI.NZ

I recommend the book by Silverman:
Density Estimation for Statistics and Data Analysis,
Chapman & Hall 1986.
************************************************************
************************************************************/

/* Input : data  - Data points for smooth
           px,py - Plotting points for smooth
   Output: estimated density */

/* Simple kernel smoothing algorithm */
proc kden_n(data,px) ;
     local h,nr,nc,sm,nx,di,n,smc,i,j ;
     nr = rows(px) ;
     nc = cols(data) ;
     sm = zeros(nr,nc) ;
     i = 1 ;
     do until(i > nc) ;
        di = sortc(packr(data[.,i]),1) ;
        n = rows(di) ;
        nx = int(coreleft./(40.*n)) ;
        h = 1./stdc(di)./(n^0.2) ;
        clear smc ;
        j = 1 ;
        do until(j > nr) ;
           smc = smc | h.*meanc(pdfn(h.*(px[j:minc((j+nx-1)|nr)]' - di))) ;
           j = j + nx ;
        endo ;
        sm[.,i] = trimr(smc,1,0) ;
        i = i + 1 ;
     endo ;
     retp(sm) ;
endp ;

/* Adaptive Kernel smoothing algorithm */
proc kden_na(data,px) ;
     local l,h,nr,nc,sm,nx,di,n,smc,i,j ;
     nr = rows(px) ;
     nc = cols(data) ;
     sm = zeros(nr,nc) ;
     i = 1 ;
     do until(i > nc) ;
        di = sortc(packr(data[.,i]),1) ;
        n  = rows(di) ;
        if n < 200 ;
           l = kden_n(di,di) ;
        else ;
           l = pdfn(di) ;
        endif ;
        l  = 1./sqrt(l./exp(meanc(ln(l)))) ;
        nx = int(coreleft./(40.*n)) ;
        h = 1.1/(l.*stdc(di)./(n^0.2)) ;
        clear smc ;
        j = 1 ;
        do until(j > nr) ;
           smc = smc | meanc(h.*pdfn(h.*(px[j:minc(j+nx-1|nr)]' - di))) ;
           j = j + nx ;
        endo ;
        sm[.,i] = trimr(smc,1,0) ;
        i = i + 1 ;
     endo ;
     retp(sm) ;
endp ;

/* 2 Dimensional kernel smooth */
proc kden_n2(data,px,py) ;
     local h,nr,nc,sm,nx,xy,x,y,n,r,c,dy,dxy,smc,i,j,k ;
     if cols(data) /= 2 ;
        errorlog " Error: KDEN_N2 - Data must have 2 columns" ;
        retp(error(1)) ;
     endif ;
     px = vecr(px) ;
     py = vecr(py)' ;
     nr = rows(px) ;
     nc = cols(py) ;
     sm = zeros(nr,nc) ;
     xy = packr(sortc(data,1)) ;
     n  = rows(xy) ;
     x = xy[.,1] ;
     y = xy[.,2] ;
     h  = (8.*n^0.2)./(2.*stdc(xy)) ;
     clear xy ;
     x  = sqrt(h[1]).*x  ;
     px = sqrt(h[1]).*px ;
     y  = sqrt(h[2]).*y  ;
     py = sqrt(h[2]).*py ;
     nx = int(coreleft./(40.*n)) ;
     print " Kernel smoothing " ;;
     r = csrlin ; c = csrcol ;
     print "  0.0 % completed" ;;
     format /rd 5,1 ;
     i = 1 ;
     do until(i > nc) ;
        clear smc ;
        dy = -(y - py[i])^2 ;
        j = 1 ;
        do until(j > nr) ;
           locate r,c ; print (100.*((i-1) + j./nr)./nc) ;;
           k   = minc(j+nx-1|nr) ;
           dxy = dy - (x - px[j:k]')^2 ;
           smc = smc | sumc(exp(dxy)) ;
           j = j + nx ;
        endo ;
        sm[i,.] = trimr(smc,1,0)' ;
        i = i + 1 ;
     endo ;
     locate r,c ; print 100 ;; print "\l" ;
     retp(((4.*prodc(h))./(n.*sqrt(2.*pi))).*sm) ;
endp ;

/* Note the smoothing parameter/bandwidth is chosen to be
   that which is optimal for a sample from a normal distribution
   You could add another input parameter to the functions which
   multiplies h to adjust the amount of smoothing manually */
