Home > wafo > kdetools > hldpi.m

hldpi

PURPOSE ^

L-stage Direct Plug-In estimate of smoothing parameter.

SYNOPSIS ^

h=hldpi(A,kernel,L,inc);

DESCRIPTION ^

 HLDPI L-stage Direct Plug-In estimate of smoothing parameter. 
  
  CALL: hs = hldpi(data,kernel,L) 
  
        hs = one dimensional value for smoothing parameter 
             given the data and kernel.  size 1 x D 
    data   = data matrix, size N x D (D = # dimensions ) 
    kernel = 'epanechnikov'  - Epanechnikov kernel. 
             'biweight'      - Bi-weight kernel. 
             'triweight'     - Tri-weight kernel. 
             'triangluar'    - Triangular kernel. 
             'gaussian'      - Gaussian kernel 
             'rectangular'   - Rectanguler kernel. 
             'laplace'       - Laplace kernel. 
             'logistic'      - Logistic kernel. 
         L = 0,1,2,3,...   (default 2) 
  
   Note that only the first 4 letters of the kernel name is needed. 
  
   Example: 
    x  = wnormrnd(0,1,50,1); 
    hs = hldpi(x,'gauss',1); 
  
  See also  hste, hbcv, hboot, hos, hlscv, hscv, hstt, kde, kdefun

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function h=hldpi(A,kernel,L,inc); 
002 %HLDPI L-stage Direct Plug-In estimate of smoothing parameter. 
003 % 
004 % CALL: hs = hldpi(data,kernel,L) 
005 % 
006 %       hs = one dimensional value for smoothing parameter 
007 %            given the data and kernel.  size 1 x D 
008 %   data   = data matrix, size N x D (D = # dimensions ) 
009 %   kernel = 'epanechnikov'  - Epanechnikov kernel. 
010 %            'biweight'      - Bi-weight kernel. 
011 %            'triweight'     - Tri-weight kernel. 
012 %            'triangluar'    - Triangular kernel. 
013 %            'gaussian'      - Gaussian kernel 
014 %            'rectangular'   - Rectanguler kernel. 
015 %            'laplace'       - Laplace kernel. 
016 %            'logistic'      - Logistic kernel. 
017 %        L = 0,1,2,3,...   (default 2) 
018 % 
019 %  Note that only the first 4 letters of the kernel name is needed. 
020 % 
021 %  Example: 
022 %   x  = wnormrnd(0,1,50,1); 
023 %   hs = hldpi(x,'gauss',1); 
024 % 
025 % See also  hste, hbcv, hboot, hos, hlscv, hscv, hstt, kde, kdefun 
026  
027 %  Wand,M.P. and Jones, M.C. (1995) 
028 % 'Kernel smoothing' 
029 %  Chapman and Hall, pp 67--74 
030  
031  
032 %tested on: matlab 5.2 
033 % revised pab Agu 2005 
034 % -fixed a bug for kernels other than Gaussian 
035 % -made it more general: L>3 is now possible 
036 % revised pab nov2004 
037 % - kernel2 = 'gaus' 
038 % -added nargchk 
039 % revised pab Dec2003 
040 % revised pab 16.10.1999 
041 %  updated to matlab 5.x changed name from hdpi -> hldpi 
042 %  changed the gridding -> hldpi is now faster 
043 % taken from kdetools             Christian C. Beardah 1994 
044  
045  
046 error(nargchk(1,4,nargin)) 
047 if nargin<2|isempty(kernel) 
048   kernel='gauss'; 
049 end 
050 if nargin<3|isempty(L), 
051   L=2; 
052 else 
053   L=abs(L); 
054 end; 
055 if nargin<4|isempty(inc) 
056   inc=128; 
057 end 
058  
059 nfft = inc*2; 
060 [n, d] = size(A); 
061 if (n==1) & (d>1), 
062   A=A.'; 
063   n=d; 
064   d=1; 
065 end 
066  
067 xmin = min(A);    % Find the minimum value of A. 
068 xmax = max(A);    % Find the maximum value of A. 
069 xrange = xmax-xmin; % Find the range of A. 
070  
071 sigmaA = std(A); 
072 iqr = abs(diff(qlevels2(A,[75 25]),1,1)); % interquartile range 
073 k = find(iqr>0); 
074 if any(k) 
075   sigmaA(k) = min(sigmaA(k), iqr(k)/1.349); 
076 end 
077  
078 % R= int(mkernel(x)^2) 
079 % mu2= int(x^2*mkernel(x)) 
080 [mu2,R] = kernelstats(kernel); 
081  
082 STEconstant = R /(mu2^(2)*n); 
083  
084 % xa holds the x 'axis' vector, defining a grid of x values where 
085 % the k.d. function will be evaluated and plotted. 
086  
087 ax1 = xmin-xrange/8; 
088 bx1 = xmax+xrange/8; 
089  
090 kernel2 = 'gaus'; 
091 [mu2,R] = kernelstats(kernel2); % Bug fix: Do not delete this: pab aug2005 
092  
093 for dim=1:d 
094   s = sigmaA(dim); 
095   ax = ax1(dim); 
096   bx = bx1(dim); 
097  
098   xa = linspace(ax,bx,inc).'; 
099   xn = linspace(0,bx-ax,inc); 
100  
101   c = gridcount(A(:,dim),xa); 
102  
103   %c=zeros(inc,1); 
104   %deltax =(bx-ax)/(inc-1); 
105   %xn     = (0:(inc-1))*deltax; 
106   %binx   = floor((A(:,dim)-ax)/deltax)+1; 
107  
108  
109   % Obtain the grid counts. 
110   %c = full(sparse(binx,1,(xa(binx+1)-A(:,dim)),inc,1)+... 
111   %         sparse(binx+1,1,(A(:,dim)-xa(binx)),inc,1))/deltax; 
112  
113   r   = 2*L+4; 
114   rd2 = L+2; 
115  
116   % Eq. 3.7 in Wand and Jones (1995) 
117   PSI_r =  (-1)^(rd2)*prod(rd2+1:r)/(sqrt(pi)*(2*s)^(r+1)); 
118   PSI   = PSI_r; 
119   if L>0 
120     % High order derivatives of the Gaussian kernel 
121     [Kd{1:L}] = deriv(0,kernel2); 
122      
123     % L-stage iterations to estimate PSI_4 
124     for ix = L:-1:1 
125       gi = (-2*Kd{ix}/(mu2*PSI*n))^(1/(2*ix+5)); 
126  
127       % Obtain the kernel weights. 
128       [KW0{1:ix}] = deriv(xn/gi,kernel2); 
129       kw=[KW0{ix},0,KW0{ix}([inc:-1:2])].'; % Apply 'fftshift' to kw. 
130  
131       % Perform the convolution. 
132       z=real(ifft(fft(c,nfft).*fft(kw))); 
133  
134       PSI = sum(c.*z(1:inc))/(n^2*gi^(2*ix+3)); 
135     end 
136   end 
137   h(dim)=(STEconstant/PSI)^(1/5); 
138 end

Mathematical Statistics
Centre for Mathematical Sciences
Lund University with Lund Institute of Technology

Comments or corrections to the WAFO group


Generated on Thu 06-Oct-2005 02:21:16 for WAFO by m2html © 2003