HLDPI2FFT L-stage DPI estimate of smoothing parameter for 2D data (DPI=Direct Plug-In) CALL h=hldpi2fft(A,kernel,L) hs = smoothing parameter data = data matrix, size N x D (D = # dimensions ==2) kernel = 'gaussian' - Gaussian kernel 'epanechnikov' - Epanechnikov kernel. 'biweight' - Bi-weight kernel. 'triweight' - Tri-weight kernel. 'triangluar' - Triangular 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,2); hs = hldpi2fft(x,'gauss',1); See also hste, hbcv, hboot, hos, hlscv, hscv, hstt, kde, kdefun
High order partial derivatives of the Gaussian kernel. | |
D-dimensional histogram using linear binning. | |
Return 2'nd order moment of kernel pdf | |
Correlation coefficients. | |
Display message and abort function. | |
Two-dimensional discrete Fourier Transform. | |
Shift zero-frequency component to center of spectrum. | |
Two-dimensional inverse discrete Fourier transform. | |
Linearly spaced vector. | |
X and Y arrays for 3-D plots. | |
Standard deviation. | |
Compare strings ignoring case. |
001 function h=hldpi2fft(A,kernel,L) 002 %HLDPI2FFT L-stage DPI estimate of smoothing parameter for 2D data 003 % (DPI=Direct Plug-In) 004 % 005 % CALL h=hldpi2fft(A,kernel,L) 006 % 007 % hs = smoothing parameter 008 % data = data matrix, size N x D (D = # dimensions ==2) 009 % kernel = 'gaussian' - Gaussian kernel 010 % 'epanechnikov' - Epanechnikov kernel. 011 % 'biweight' - Bi-weight kernel. 012 % 'triweight' - Tri-weight kernel. 013 % 'triangluar' - Triangular 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,2); 023 % hs = hldpi2fft(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. (1994) 028 % "Multivariate Plug-In Bandwidth Selection", 029 % Computational statistics 030 031 % tested on: matlab 5.2 032 % history: 033 % revised pab aug2005 034 % -added ad hoc support for other kernels than Gaussian 035 % revised pab dec2003 036 % -fixed some bugs 037 % -added todo comments 038 % -added call to gridcount 039 % -updated example 040 % -fixed a bug 041 % revised pab 3.11.1999 042 % from kdetools by Christian C. Beardah 1996 043 044 045 046 047 [n,d]=size(A); 048 049 if d~=2, 050 error('Wrong size of data') 051 end 052 if nargin<2 |isempty(kernel) 053 kernel = 'gauss'; 054 end 055 056 if nargin<3|isempty(L), 057 L=2; 058 else 059 L=min(abs(L),3); 060 end; 061 062 if any(strcmpi(kernel(1:4),{'gaus','norm'})), 063 r = 4; 064 else, 065 r=1; 066 end; 067 068 inc=128; 069 070 P=2*inc; 071 072 xmin = min(A,[],1); % Find the minimum value of x. 073 xmax = max(A,[],1); % Find the maximum value of x. 074 xrange=xmax-xmin; % Find the range of x. 075 076 % xa holds the x 'axis' vector, defining a grid of x values where 077 % the k.d. function will be evaluated and plotted. 078 079 ax=xmin-xrange/4; 080 bx=xmax+xrange/4; 081 082 deltax = (bx-ax)/(inc-1); 083 084 X =[linspace(ax(1),bx(1),inc).' , linspace(ax(2),bx(2),inc).']; 085 086 % Find the binned kernel weights, c. 087 c = gridcount(A,X); 088 089 % R= int(mkernel(x)^2) 090 % mu2= int(x^2*mkernel(x)) 091 kernel2 = 'gauss'; 092 % values for gaussian kernel 093 [mu2ns,Rns] = kernelstats(kernel2); 094 Rns = Rns^d; 095 STEconstant2 = Rns /(mu2ns^(2)*n); 096 097 [mu2,R] = kernelstats(kernel); 098 R = R^d; 099 STEconstant = R /(mu2^(2)*n); 100 101 102 103 104 K22=1/(2*pi); 105 K40=3/(2*pi); 106 K04=3/(2*pi); 107 108 K42=-3/(2*pi); 109 K24=-3/(2*pi); 110 K60=-15/(2*pi); 111 K06=-15/(2*pi); 112 113 K44=9/(2*pi); 114 K62=15/(2*pi); 115 K26=15/(2*pi); 116 K80=105/(2*pi); 117 K08=105/(2*pi); 118 119 sigma1=std(A(:,1)); 120 sigma2=std(A(:,2)); 121 122 rho = corrcoef(A); 123 rho = rho(1,2); 124 125 if L==3, 126 127 lam100=945; 128 lam010=945; 129 lam82=105+840*rho^2; 130 lam28=105+840*rho^2; 131 lam64=45+540*rho^2+360*rho^4; 132 lam46=45+540*rho^2+360*rho^4; 133 134 psi100=lam100/(128*pi*sigma1^11*sigma2*(1-rho^2)^5.5); 135 psi010=lam010/(128*pi*sigma1*sigma2^11*(1-rho^2)^5.5); 136 psi82=lam82/(128*pi*sigma1^9*sigma2^3*(1-rho^2)^5.5); 137 psi28=lam28/(128*pi*sigma1^3*sigma2^9*(1-rho^2)^5.5); 138 psi64=lam64/(128*pi*sigma1^7*sigma2^5*(1-rho^2)^5.5); 139 psi46=lam46/(128*pi*sigma1^5*sigma2^7*(1-rho^2)^5.5); 140 141 % The following are the smoothing parameters for 142 % estimation of psi80, psi08, psi62, psi26 and psi44: 143 144 a80=(-2*K80/(mu2ns*(psi100+psi82)*n))^(1/12); 145 a08=(-2*K08/(mu2ns*(psi28+psi010)*n))^(1/12); 146 a62=(-2*K62/(mu2ns*(psi82+psi64)*n))^(1/12); 147 a26=(-2*K26/(mu2ns*(psi46+psi28)*n))^(1/12); 148 a44=(-2*K44/(mu2ns*(psi64+psi46)*n))^(1/12); 149 150 psi80=psim('80',A,a80); 151 psi08=psim('08',A,a08); 152 psi62=psim('62',A,a62); 153 psi26=psim('26',A,a26); 154 psi44=psim('44',A,a44); 155 156 end; 157 158 if L==2 | L==3, 159 160 if L==2, 161 lam80=105; 162 lam08=105; 163 lam62=15+90*rho^2; 164 lam26=15+90*rho^2; 165 lam44=9+72*rho^2+24*rho^4; 166 167 psi80=lam80/(64*pi*sigma1^9*sigma2*(1-rho^2)^4.5); 168 psi08=lam08/(64*pi*sigma1*sigma2^9*(1-rho^2)^4.5); 169 psi62=lam62/(64*pi*sigma1^7*sigma2^3*(1-rho^2)^4.5); 170 psi26=lam26/(64*pi*sigma1^3*sigma2^7*(1-rho^2)^4.5); 171 psi44=lam44/(64*pi*sigma1^5*sigma2^5*(1-rho^2)^4.5); 172 end; 173 174 % The following are the smoothing parameters for 175 % estimation of psi60, psi06, psi42 and psi24: 176 177 a60=(-2*K60/(mu2ns*(psi80+psi62)*n))^0.1; 178 a06=(-2*K06/(mu2ns*(psi26+psi08)*n))^0.1; 179 a42=(-2*K42/(mu2ns*(psi62+psi44)*n))^0.1; 180 a24=(-2*K24/(mu2ns*(psi44+psi26)*n))^0.1; 181 182 psi60=psim('60',A,a60); 183 psi06=psim('06',A,a06); 184 psi42=psim('42',A,a42); 185 psi24=psim('24',A,a24); 186 187 end; 188 189 if L==1 | L==2 | L==3, 190 191 if L==1, 192 lam06=15; 193 lam60=15; 194 lam24=3+12*rho^2; 195 lam42=3+12*rho^2; 196 197 psi60=-lam60/(32*pi*sigma1^7*sigma2*(1-rho^2)^3.5); 198 psi06=-lam06/(32*pi*sigma1*sigma2^7*(1-rho^2)^3.5); 199 psi42=-lam42/(32*pi*sigma1^5*sigma2^3*(1-rho^2)^3.5); 200 psi24=-lam24/(32*pi*sigma1^3*sigma2^5*(1-rho^2)^3.5); 201 end; 202 203 % The following are the smoothing parameters for 204 % estimation of psi40, psi04 and psi22: 205 206 a40=(-2*K40/(mu2ns*(psi60+psi42)*n))^0.125; 207 a04=(-2*K04/(mu2ns*(psi24+psi06)*n))^0.125; 208 a22=(-2*K22/(mu2ns*(psi42+psi24)*n))^0.125; 209 210 % Brute force technique: 211 212 psi40=psim('40',A,a40); 213 psi04=psim('04',A,a04); 214 psi22=psim('22',A,a22); 215 216 % Fourier technique: 217 218 % Psi40: 219 220 h=real(a40); 221 222 % Obtain the kernel weights. 223 224 L1=max(floor(r*h./deltax)); 225 226 L2=min(L1,inc-1); 227 228 kw=zeros(2*inc,2*inc); 229 230 [X1,Y1]=meshgrid(deltax(1)*[-L2:L2]/h,deltax(2)*[-L2:L2]/h); 231 idx(1:d) = {(inc-L2+1):(inc+L2+1)}; 232 233 kw(idx{:}) = deriv2(X1,Y1,'40')/(n*h^2); 234 235 % Apply 'fftshift' to kw. 236 237 kw = fftshift(kw); 238 239 % Perform the convolution. 240 241 z = real(ifft2(fft2(c,P,P).*fft2(kw))); 242 243 z=z(1:inc,1:inc).*(z(1:inc,1:inc)>0); 244 245 psi40=sum(sum(c.*z(1:inc,1:inc)))/(n^2*h^5); 246 247 end; 248 249 if L==0, 250 251 lam04=3; 252 lam40=3; 253 lam22=3/2; 254 255 psi40=-lam40/(16*pi*sigma1^5*sigma2*(1-rho^2)^2.5); 256 psi04=-lam04/(16*pi*sigma1*sigma2^5*(1-rho^2)^2.5); 257 psi22=-lam22/(16*pi*sigma1^3*sigma2^3*(1-rho^2)^2.5); 258 259 end; 260 261 262 263 h1=(psi04^0.75*R/(n*mu2^2*psi40^0.75*(sqrt(psi40*psi04)+psi22)))^(1/6); 264 h2=h1*(psi40/psi04)^0.25; 265 266 h=real([h1, h2]); 267 268 269 return 270 271 272 function p=psim(m,A,h) 273 274 % PSIM Calculate psi_m by brute force. 275 % 276 % PSIM(M,A,H) 277 % 278 % Christian C. Beardah 1996 279 280 n=length(A); 281 282 Xi=zeros(n,n); 283 Xj=Xi; 284 Yi=Xi; 285 Yj=Xi; 286 287 o=ones(1,n); 288 289 Xj=A(:,1)*o; 290 Xi=Xj'; 291 Yj=A(:,2)*o; 292 Yi=Yj'; 293 294 p=sum(sum(deriv2((Xi-Xj)/h,(Yi-Yj)/h,m) ))/(h*n^2); 295 return 296
Comments or corrections to the WAFO group