function hout = bwdipi(data,ihtype) 
% BWDIPI, Devroye's Improved Plug In, for 1-d kernel density estimation
%     As described in Sections 9 & 10 of Devroye's 1997 TEST discussion
%     paper, attributed to Berlinet and Devroye 1994, (Pub. Inst. 
%                   Stat. Univ. Paris)
% Inputs:
%     data   - column vector of data,
%     ihtype - 0 all 3 types
%            - 1 Improved Plug In (default, Devroye's recco)
%            - 2 Double Kernel
%            - 3 L1 Plug In
% Output:
%     h - Devroye's Improved Plug In bandwidth
%

%    Copyright (c) J. S. Marron 1997


%  Set parameters and defaults according to number of input arguments
%
if nargin == 1 ;    %  only 1 argument input, use default ihtype
  iihtype = 1 ;
else ;              %  xgrid was specified, use that
  iihtype = ihtype ;
end ;



%  First calculate common quantities
%
ndat = length(data) ;
mindat = min(data) ;
maxdat = max(data) ;
ngrid = 200 ;      %  Grid size for numerical integrals

sighat3 = std(data) ;
          %  using formula (3) from Devroye's paper, based on sample sd

sdata = sort(data) ;
x1q = sdata(floor(ndat / 4)) ;
x3q = sdata(floor(3 * ndat / 4)) ;
sighat4 = (x3q - x1q) / 1.34898 ;
          %  using formula (4) from Devroye's paper, based on IQR



%  Calculate h L1 Plug-In, if requested
%
if iihtype == 3  |  iihtype == 0 ;

  %  Each block is one step in Devroye's paper

  hp = 2.279 * sighat4 * ndat^(-1/5) ; 
  %  this the "L1 reference bandwidth", href,L1

    left = mindat - 1.5 * hp ;    %  use larger grid since 
                                  %  L is supported [-1.5,1.5]
    right = maxdat + 1.5 * hp ;
    del = (right - left) / (ngrid - 1) ;
    xgrid = linspace(left,right,ngrid)' ;
    arg = (vec2mat(xgrid',ndat) - vec2mat(data,ngrid)) / hp ;
    fnh = sum((abs(arg) < 1) .* (1 - arg.^2)) ;
    fnh = (3/4) * fnh' / (ndat * hp) ;
  a = del * sum(sqrt(fnh)) ;

    arg = arg / 1.5 ;
    gnh = (abs(arg) < .5) .* (7 - 31 * arg.^2) ;
    gnh = gnh + (abs(arg) >= .5) .* (abs(arg) < 1) .* (arg.^2 - 1) ;
    gnh = sum(gnh) ;
    gnh = (1/4) * gnh' / (ndat * 1.5 * hp) ;
    intfmg = del * sum(abs(fnh - gnh)) ;
    rintkml2 = 935 / 2247 ;
            %  sqrt(int((K-L)^2)), from devker.m
  r = a * rintkml2 / (sqrt(ndat * hp) * intfmg);

  hpp = hp * max([1; (10 * r)^(2/5)]) ;
  
    left = mindat - 1.5 * hpp ;    %  use larger grid since 
                                  %  L is supported [-1.5,1.5]
    right = maxdat + 1.5 * hpp ;
    del = (right - left) / (ngrid - 1) ;
    xgrid = linspace(left,right,ngrid)' ;
    arg = (vec2mat(xgrid',ndat) - vec2mat(data,ngrid)) / hpp ;
    fnh = sum((abs(arg) < 1) .* (1 - arg.^2)) ;
    fnh = (3/4) * fnh' / (ndat * hpp) ;
    arg = arg / 1.5 ;
    gnh = (abs(arg) < .5) .* (7 - 31 * arg.^2) ;
    gnh = gnh + (abs(arg) >= .5) .* (abs(arg) < 1) .* (arg.^2 - 1) ;
    gnh = sum(gnh) ;
    gnh = (1/4) * gnh' / (ndat * 1.5 * hpp) ;
    intfmg = del * sum(abs(fnh - gnh)) ;
    intx2k = 1 / 5 ;
  b = 2 * intfmg / (hpp^2 * intx2k) ;

  hmsL1 = 2.71042 * sighat3 * ndat^(-1/5) ;

  minarg1 = (sqrt(15 / (2 * pi)) * a / b)^(2/5) * ndat^(-1/5) ;
  hl1pi = min([minarg1; hmsL1]) ;

end ;



%  Next calculate h Double kernel, if needed
%
if iihtype <= 2  ;

  hrefl1 = 1.6644 * sighat4 * ndat^(-1/5) ; 
  %  this the "L1 reference bandwidth", href,l1

    left = mindat - 1.5 * hrefl1 ;    %  use larger grid since 
                                  %  L is supported [-1.5,1.5]
    right = maxdat + 1.5 * hrefl1 ;
    del = (right - left) / (ngrid - 1) ;
    xgrid = linspace(left,right,ngrid)' ;
    arg = (vec2mat(xgrid',ndat) - vec2mat(data,ngrid)) / hrefl1 ;
    fnh = sum((abs(arg) < 1) .* (1 - arg.^2)) ;
    fnh = (3/4) * fnh' / (ndat * hrefl1) ;
    arg = arg / 1.5 ;
    gnh = (abs(arg) < .5) .* (7 - 31 * arg.^2) ;
    gnh = gnh + (abs(arg) >= .5) .* (abs(arg) < 1) .* (arg.^2 - 1) ;
    gnh = sum(gnh) ;
    gnh = (1/4) * gnh' / (ndat * 1.5 * hrefl1) ;
  jnhrefl1 = del * sum(abs(fnh - gnh)) ;

  %  Find limit a
  %
  intkml = 447 / 775 ;
          %  int(|K - L|), calculated numerically by devker.m
  a = hrefl1 ;
    flag = (sdata(1:(ndat-2)) + 2 * a <= sdata(2:(ndat-1))) .* ...
                     (sdata(2:(ndat-1)) <= sdata(3:ndat) - 2 * a) ;
    flag = [(sdata(1) + 2 * a <= sdata(2)) ; flag] ;
    flag = [flag; ((sdata(ndat-1) + 2 * a) <= sdata(ndat))] ;
  chia = (intkml / ndat) * sum(flag) ;
  while chia < jnhrefl1
    a = a / 2 ;
      flag = (sdata(1:(ndat-2)) + 2 * a <= sdata(2:(ndat-1))) .* ...
                     (sdata(2:(ndat-1)) <= sdata(3:ndat) - 2 * a) ;
      flag = [(sdata(1) + 2 * a <= sdata(2)) ; flag] ;
      flag = [flag; ((sdata(ndat-1) + 2 * a) <= sdata(ndat))] ;
    chia = (intkml / ndat) * sum(flag) ;
  end ;

  %  Find limit b
  %
  c = 864 / 373 ;
          %  Lipschitz constant for K - L, calculated by devker.m
  b = hrefl1 ;
  muhat = mean(data) ;
    arg = abs(data - muhat) ;
    arg = (2 * b + arg) .* c .* arg ./ b^2;
    arg = min([arg'; (2 * intkml * ones(1,ndat))])' ;
  xib = intkml - (1 / ndat) * sum(arg) ;
  while xib < jnhrefl1
    b = b * 2 ;
      arg = abs(data - muhat) ;
      arg = (2 * b + arg) .* c .* arg ./ b^2;
      arg = min([arg'; (2 * intkml * ones(1,ndat))])' ;
    xib = intkml - (1 / ndat) * sum(arg) ;
  end ;
  
    hacc = 1.2 ;
          %  need to carefully think about this.  Ask Luc?
    nhgrid = log(b/a) / log(hacc) + 1 ;
  hgrid = logspace(log10(a),log10(b),nhgrid) ;
  vjnh = [] ;
  for i = 1:nhgrid;
    hi = hgrid(i) ;
      left = mindat - 1.5 * hi ;    %  use larger grid since 
                                        %  L is supported [-1.5,1.5]
      right = maxdat + 1.5 * hi ;
      del = (right - left) / (ngrid - 1) ;
      xgrid = linspace(left,right,ngrid)' ;
      arg = (vec2mat(xgrid',ndat) - vec2mat(data,ngrid)) / hi ;
      fnh = sum((abs(arg) < 1) .* (1 - arg.^2)) ;
      fnh = (3/4) * fnh' / (ndat * hi) ;
      arg = arg / 1.5 ;
      gnh = (abs(arg) < .5) .* (7 - 31 * arg.^2) ;
      gnh = gnh + (abs(arg) >= .5) .* (abs(arg) < 1) .* (arg.^2 - 1) ;
      gnh = sum(gnh) ;
      gnh = (1/4) * gnh' / (ndat * 1.5 * hi) ;
    jnh = del * sum(abs(fnh - gnh)) ;
    vjnh = [vjnh; jnh] ;
  end ;
  [temp,ihdk] = min(vjnh) ;
  hdk = hgrid(ihdk) ;

end ;



%  Now "improve" if requested
%
if iihtype <= 1  ;

  %  Each block is one step in Devroye's paper

  hp = hdk ; 
  %  Now start with the "double kernel bandwidth"

    left = mindat - 1.5 * hp ;    %  use larger grid since 
                                  %  L is supported [-1.5,1.5]
    right = maxdat + 1.5 * hp ;
    del = (right - left) / (ngrid - 1) ;
    xgrid = linspace(left,right,ngrid)' ;
    arg = (vec2mat(xgrid',ndat) - vec2mat(data,ngrid)) / hp ;
    fnh = sum((abs(arg) < 1) .* (1 - arg.^2)) ;
    fnh = (3/4) * fnh' / (ndat * hp) ;
  a = del * sum(sqrt(fnh)) ;

    arg = arg / 1.5 ;
    gnh = (abs(arg) < .5) .* (7 - 31 * arg.^2) ;
    gnh = gnh + (abs(arg) >= .5) .* (abs(arg) < 1) .* (arg.^2 - 1) ;
    gnh = sum(gnh) ;
    gnh = (1/4) * gnh' / (ndat * 1.5 * hp) ;
    intfmg = del * sum(abs(fnh - gnh)) ;
    rintkml2 = 935 / 2247 ;
            %  sqrt(int((K-L)^2)), from devker.m
  r = a * rintkml2 / (sqrt(ndat * hp) * intfmg);

  hpp = hp * max([1; (10 * r)^(2/5)]) ;
  
    left = mindat - 1.5 * hpp ;    %  use larger grid since 
                                  %  L is supported [-1.5,1.5]
    right = maxdat + 1.5 * hpp ;
    del = (right - left) / (ngrid - 1) ;
    xgrid = linspace(left,right,ngrid)' ;
    arg = (vec2mat(xgrid',ndat) - vec2mat(data,ngrid)) / hpp ;
    fnh = sum((abs(arg) < 1) .* (1 - arg.^2)) ;
    fnh = (3/4) * fnh' / (ndat * hpp) ;
    arg = arg / 1.5 ;
    gnh = (abs(arg) < .5) .* (7 - 31 * arg.^2) ;
    gnh = gnh + (abs(arg) >= .5) .* (abs(arg) < 1) .* (arg.^2 - 1) ;
    gnh = sum(gnh) ;
    gnh = (1/4) * gnh' / (ndat * 1.5 * hpp) ;
    intfmg = del * sum(abs(fnh - gnh)) ;
    intx2k = 1 / 5 ;
  b = 2 * intfmg / (hpp^2 * intx2k) ;

  hmsL1 = 2.71042 * sighat3 * ndat^(-1/5) ;

  minarg1 = (sqrt(15 / (2 * pi)) * a / b)^(2/5) * ndat^(-1/5) ;
  hipi = min([minarg1; hmsL1]) ;

end ;



if iihtype == 1 ;
  hout = hipi ;
elseif iihtype == 2 ;
  hout = hdk ;
elseif iihtype == 3 ;
  hout = hl1pi ;
elseif iihtype == 0 ;
  hout = [hipi; hdk; hl1pi] ;
end ;
