disp('Running MATLAB script file bwdipit.m') ;
%
%    FOR DEVELOPMENT AND TESTING OF MATLAB FUNCTION bwdipi.m,
%    Devroye's Improved Plug In BandWidth

itest = 2 ;     %  1 - initial development
                %  2 - construction of pic, and data fro message for luc


%  Generate .5 N(-3,1) + .5 N(3,1) density
%
ng = 200 ;
xg = linspace(-6,6,ng)' ;
f = .5 * (1 / sqrt(2 * pi)) * ...
               (exp(-(xg - 3).^2 / 2) + exp(-(xg + 3).^2 / 2)) ;


%  Generate .5 N(-3,1) + .5 N(3,1) data
%
n = 100 ;
rand('seed',93274923) ;
randn('seed',38964987) ;
data = randn(n,1) + 3 * (2 * (rand(n,1) > .5) - 1) ;


if itest == 1 ;

  plot(xg,f) ;
            %  put this down, then can plot diagnostic plots on top

  tic ;
  hout = bwdipi(data,1) ;
            %  recall 
            %     ihtype - 0 all 3 types
            %            - 1 Improved Plug In (default, Devroye's recco)
            %            - 2 Double Kernel
            %            - 3 L1 Plug In
  toc ;

  hout



  %  Add estimate to plot
  mindat = min(data) ;
  maxdat = max(data) ;
  left = mindat - 1.5 * hout ;    %  use larger grid since 
                                  %  L is supported [-1.5,1.5]
  right = maxdat + 1.5 * hout ;
  del = (right - left) / (ng - 1) ;
  xg = linspace(left,right,ng)' ;
  arg = (vec2mat(xg',n) - vec2mat(data,ng)) / hout ;
  fnh = sum((abs(arg) < 1) .* (1 - arg.^2)) ;
  fnh = (3/4) * fnh' / (n * hout) ;
  hold on ;
    plot(xg,fnh,'m') ;
    text(-2,.15,['h = ' num2str(hout)]) ;
  hold off ;

elseif itest == 2 ;

  %  get data driven bandwidths
  hout = bwdipi(data,0) ;
            %     ihtype - 0 all 3 types

  %  do plots
  for ip = 1:4 ;
    subplot(2,2,ip) ;

    if ip <= 3 ;
      arg = (vec2mat(xg',n) - vec2mat(data,ng)) / hout(ip) ;
      fnh = sum((abs(arg) < 1) .* (1 - arg.^2)) ;
      fnh = (3/4) * fnh' / (n * hout(ip)) ;
    else ;
      hsjpi = bwsjpib(data,[-6; 6; ng]) ;
      fnh = gpkde(data,hsjpi,[-6; 6; ng]) ;
    end ;

    if ip == 1 ;
      titstr = 'Epa. kernel, Improved Plugin' ;
      textstr = ['h = ' num2str(hout(ip))] ;
    elseif ip == 2 ;
      titstr = 'Epa. kernel, Double Kernel' ;
      textstr = ['h = ' num2str(hout(ip))] ;
    elseif ip == 3 ;              
      titstr = 'Epa. kernel, L1 Plugin' ;
      textstr = ['h = ' num2str(hout(ip))] ;
    elseif ip == 4 ;
      titstr = 'Normal kernel, SJPI' ;
      textstr = ['h = ' num2str(hsjpi)] ;
    end ;

    plot(xg,[f,fnh]) ;
      title(titstr) ;
      text(-2,.15,textstr) ;

  end ;



  %  WRITE ANSWERS INTO FILE
  fid = fopen('/matlab/stevefun/tests/bwdipit.out','wt') ;

    titlstr1 = 'Output from the MATLAB Script bwdipit.m,   ' ;
  titlstr1 = [titlstr1,date] ;
  cntbytes = fprintf(fid,'%1s\n\n\n\n',titlstr1) ;


    titlstr3 = 'Data driven bandwidths were: ' ;
  cntbytes = fprintf(fid,'%1s\n',titlstr3) ;

    forstr = '    Improved Plug In:  %12.8f\n' ;
  cntbytes = fprintf(fid,forstr,hout(1)) ;

    forstr = '       Double Kernel:  %12.8f\n' ;
  cntbytes = fprintf(fid,forstr,hout(2)) ;

    forstr = '          L1 Plug In:  %12.8f\n\n\n' ;
  cntbytes = fprintf(fid,forstr,hout(3)) ;


    titlstr3 = 'Based on the data: ' ;
  cntbytes = fprintf(fid,'%1s\n',titlstr3) ;

    forstr = '  %12.8f\n' ;
  cntbytes = fprintf(fid,forstr,data) ;


  fclose('all') ;


end ;
