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);

## CROSS-REFERENCE INFORMATION

This function calls:
 deriv 4th, 6th, 8th and 10th derivatives of the kernel function. gridcount D-dimensional histogram using linear binning. kernelstats Return 2'nd order moment of kernel pdf qlevels2 Calculates quantile levels which encloses P% of data diff Difference and approximate derivative. error Display message and abort function. fft Discrete Fourier transform. ifft Inverse discrete Fourier transform. linspace Linearly spaced vector. std Standard deviation.
This function is called by:
 kdeplot Computes and plots a kernel density estimate of PDF

## 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'
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