function [kde,xgrid,mker] = CHkde(data,vh,itype,vxgrid,eptflag) 
% CHKDE, Censored Hazard Kernel Density Estimate (1-d, Gaussian Kernel)
%     Does 1-d kernel density estimation, using a binned 
%     implementations, with the bandwidth user specified 
%     (can be vector)
%   Can use first 1, 2, 3, 4 or 5 arguments.
% Inputs:
%     data   - either n x 1 column vector of uncensored data
%                  or n x 2 matrix of censored data, with:
%                      X's in first column,
%                      delta's in second column, with values:
%                          1 - when X is the actual value
%                          0 - when X is the (right) censoring time
%                                (i.e. the actual value is only
%                                      known to be larger)
%     vh     - vector of bandwidths (required)
%    itype   - index of estimation type:
%                  1  density estimation (default)
%                  2  hazard rate estimation
%     vxgrid - vector of parameters
%                  0 (or not specified)  -  use endpts of data and 401 bins
%                  [le; lr]  -  le is left end, re is right, 401 bins
%                         (get error message and no return if le > lr)
%                  [le; lr; nb] - le left, re right, and nb bins
%    eptflag - endpoint truncation flag (only has effect when imptyp = 0):
%                  0  -  move data outside range to nearest endpoint
%                  1 (or not specified)  -  truncate data outside range
%                        CAUTION:  this default is the opposite
%                            of gplbinr, but makes sense, since
%                            often want to eliminate wierd edge
%                            effects in these settings
% Output:
%     (none)  -  Draws a graph of the result (in the current axes)
%     kde     -  col vector of heights of kernel kernel density estimate,
%                    unless vh is a vector (then have matrix, with
%                    corresponding cols as density estimates)
%     xgrid   -  col vector grid of points at which estimate(s) are 
%                    evaluted (useful for plotting), unless grid is input,
%                    can also get this from linspace(le,re,nb)'  
%     mker    -  matrix (vector) of kernel functions, evaluated at xgrid,
%                    which can be plotted to show "effective window 
%                    sizes" (currently scaled to have mass 0.05, may
%                    need some vertical rescaling)
%
% Assumes path can find personal functions:
%    vec2mat.m
%    gplbinr.m
%    CHlbinr.m
%    bwsjpib.m
%    bwos.m
%    bwrot.m
%    bwsnr.m

%    Copyright (c) J. S. Marron 2000



%  Set parameters and defaults according to number of input arguments
%
if nargin <= 2 ;    %   at most 2 arguments input, use default itype
  iitype = 1 ;
else ;              %  bandwidth was specified, use that
  iitype = itype ;
end ;

if nargin <= 3 ;    %  at most 3 arguments input, use default xgrid
  ivxgrid = 0 ;
else ;              %  xgrid was specified, use that
  ivxgrid = vxgrid ;
end ;

if nargin <= 4 ;    %  at most 4 arguments input, use default endpt trunc
  ieptflag = 1 ;    %  Default
else ;
  ieptflag = eptflag ;    %  Endpt trunc was specified, so use it
end ;



%  detect whether censored or uncensored data
%
if size(data,2) == 1 ;    %  Then is uncensored
  xdat = data(:,1) ;
  vdel = ones(length(xdat),1) ;
  icen = 0 ;
else ;                    %  Then assume censored
  xdat = data(:,1) ;
  vdel = data(:,2) ;
  icen = 1 ;
end ;
n = length(xdat) ;



%  Calculate curve estimate
%
%  Specify grid parameters
nbin = 401 ;         %  Default
lend = min(xdat) ;   %  Default
rend = max(xdat) ;   %  Default
if (length(ivxgrid) == 2) | (length(ivxgrid) == 3) ;
                                 %  then use input endpoints
  lend = ivxgrid(1) ;
  rend = ivxgrid(2) ;
end ;
if length(ivxgrid) == 3 ;          %  then use number of grid points
  nbin = ivxgrid(3) ;
end ;


if lend > rend ;    %  Then bad range has been input
  disp('!!!   Error in CHkde: invalid range chosen  !!!') ;
  bincts = [] ;
else ;

  if iitype == 1  &  icen == 0  ;  
    bincts = gplbinr(data,[lend,rend,nbin],ieptflag) ;
  else ;

    [svdata, vsortind] = sort(xdat) ;
    svdel = vdel(vsortind) ;

    if icen == 1 ;    %  then are doing censored estimation
      vg = KMcdf(svdata,1 - svdel,1) ;
      vgbar = 1 - vg  ;  
    else ;
      vgbar = ones(n,1) ;
    end ;

    if iitype == 2 ;     % then are doing hazard est.
      vf = KMcdf(svdata,svdel,1) ;
      vfbar = 1 - vf  ;  
    else ;
      vfbar = ones(n,1) ;
    end ;

    vlbar = vfbar .* vgbar ;

    data = [svdata, svdel, vlbar] ;

    bincts = CHlbinr(data,[lend,rend,nbin],ieptflag) ;

  end ;

end ;


if ieptflag == 1 ;
  n = sum(xdat >= lend  &  xdat <= rend) ;
          %  put this here in case of truncations during binning
end ;


if n == 0 ;    %  Possible when all data points are censored

  kde = zeros(nbin,length(vh)) ;

else ;

  %  Loop through bandwidths
  kde = [] ;
  for ih = 1:length(vh) ;
    h = vh(ih) ;

    %  Create vector of kernel values, at equally spaced grid
    delta = (rend - lend) / (nbin - 1) ;    %  binwidth
    k = nbin - 1 ;    %  index of last nonzero entry of kernel vector
    arg = linspace(0,k * delta / h,k + 1)' ;
    kvec = exp(-(arg.^2) / 2) / sqrt(2 * pi) ;
    kvec = [flipud(kvec(2:k+1)); kvec] ;

    %  Do actual kernel density estimation
    kdeh = conv(bincts,kvec) ;
    kdeh = kdeh(k+1:k+nbin) / (n * h) ;



    kde = [kde kdeh] ;
  end ;

end ;


xgrid = linspace(lend,rend,nbin)' ;




%  Create matrix of kernels, if this is needed
%
if nargout == 3 ;
  cent = mean([lend; rend]) ;
          %  centerpoint of evaluation grid
  if length(vh) > 1 ;
    mih = vec2mat(1 ./ vh',nbin) ;
    mker = vec2mat((xgrid - cent),length(vh)) .* mih;
          %  argument of gaussian kernel
  else ;
    mih = 1 / vh ;
    mker = (xgrid - cent) .* mih;
          %  argument of gaussian kernel
  end ;
  mker = exp(-mker.^2 / 2) .* mih / sqrt(2 * pi) ;
          %  Gaussian kernels with mass 1
  mker = 0.05 * mker ;
          %  Make masses = 0.05
end ;



%  Make plots if no numerical output requested
%
if nargout == 0 ;  %  Then no numerical output, but make a plot
  plot(xgrid,kde) ;
end ;

