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
4th, 6th, 8th and 10th derivatives of the kernel function. | |
D-dimensional histogram using linear binning. | |
Return 2'nd order moment of kernel pdf | |
Calculates quantile levels which encloses P% of data | |
Difference and approximate derivative. | |
Display message and abort function. | |
Discrete Fourier transform. | |
Inverse discrete Fourier transform. | |
Linearly spaced vector. | |
Standard deviation. |
Computes and plots a kernel density estimate of PDF |
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
Comments or corrections to the WAFO group