Home > wafo > kdetools > hldpi2fft.m

hldpi2fft

PURPOSE ^

L-stage DPI estimate of smoothing parameter for 2D data

SYNOPSIS ^

h=hldpi2fft(A,kernel,L)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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

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