Home > wafo > kdetools > hldpi2.m

hldpi2

PURPOSE ^

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

SYNOPSIS ^

h=hldpi2(A,kernel,L)

DESCRIPTION ^

 HLDPI2  L-stage DPI estimate of smoothing parameter for 2D data 
          (DPI=Direct Plug-In) 
  
  CALL: hs = hldpi2(data,kernel,L) 
  
        hs = smoothing parameter 
    data   = data matrix, size N x D (D = # dimensions ==2) 
    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,2); 
    hs = hldpi2(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=hldpi2(A,kernel,L) 
002 %HLDPI2  L-stage DPI estimate of smoothing parameter for 2D data 
003 %         (DPI=Direct Plug-In) 
004 % 
005 % CALL: hs = hldpi2(data,kernel,L) 
006 % 
007 %       hs = smoothing parameter 
008 %   data   = data matrix, size N x D (D = # dimensions ==2) 
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,2); 
023 %   hs = hldpi2(x,'gauss',1); 
024 % 
025 % See also  hste, hbcv, hboot, hos, hlscv, hscv, hstt, kde, kdefun 
026  
027 % References   
028 %  Wand, M.P. and Jones, M.C. (1994)  
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 Nov2004   
036 %  added nargchk   
037 % revised pab dec2003   
038 %  -added psim   
039 % revised pab 3.11.1999 
040 %  renamed from h2dpi to hldpi2 
041 % from kdetools              Christian C. Beardah 1996 
042  
043  
044 error(nargchk(1,3,nargin))   
045 [n, d]=size(A); 
046 if d~=2, 
047   error('Wrong size of data') 
048 end 
049 if nargin<2 |isempty(kernel) 
050   kernel = 'gauss'; 
051 end 
052 if nargin<3|isempty(L), 
053   L=2; 
054 else 
055   L=min(abs(L),3); 
056 end; 
057  
058 % Transform A so rows have st. dev.=1: 
059 oldA=A; 
060 A=(oldA-ones(n,1)*mean(oldA))*diag([1./std(oldA)]); % centre data about 
061  
062 % its mean 
063  
064 K22=1/(2*pi); 
065 K40=3/(2*pi); 
066 K04=3/(2*pi); 
067  
068 K42=-3/(2*pi); 
069 K24=-3/(2*pi); 
070 K60=-15/(2*pi); 
071 K06=-15/(2*pi); 
072  
073 K44=9/(2*pi); 
074 K62=15/(2*pi); 
075 K26=15/(2*pi); 
076 K80=105/(2*pi); 
077 K08=105/(2*pi); 
078  
079 % R= int(mkernel(x)^2) 
080 % mu2= int(x^2*mkernel(x)) 
081 kernel2 = 'gauss'; 
082 % values for gaussian kernel 
083 [mu2ns,Rns] = kernelstats(kernel2); 
084 Rns = Rns^d; 
085 STEconstant2 = Rns /(mu2ns^(2)*n); 
086  
087 [mu2,R] = kernelstats(kernel); 
088 R = R^d; 
089  
090 STEconstant = R /(mu2^(2)*n); 
091  
092 sigma1=std(A(:,1)); 
093 sigma2=std(A(:,2)); 
094  
095 rho=corrcoef(A); 
096 rho=rho(1,2); 
097  
098 if L==3, 
099  
100   lam100=945; 
101   lam010=945; 
102   lam82=105+840*rho^2; 
103   lam28=105+840*rho^2; 
104   lam64=45+540*rho^2+360*rho^4; 
105   lam46=45+540*rho^2+360*rho^4; 
106  
107   psi100=lam100/(128*pi*sigma1^11*sigma2*(1-rho^2)^5.5); 
108   psi010=lam010/(128*pi*sigma1*sigma2^11*(1-rho^2)^5.5); 
109   psi82=lam82/(128*pi*sigma1^9*sigma2^3*(1-rho^2)^5.5); 
110   psi28=lam28/(128*pi*sigma1^3*sigma2^9*(1-rho^2)^5.5); 
111   psi64=lam64/(128*pi*sigma1^7*sigma2^5*(1-rho^2)^5.5); 
112   psi46=lam46/(128*pi*sigma1^5*sigma2^7*(1-rho^2)^5.5); 
113  
114 % The following are the smoothing parameters for  
115 % estimation of psi80, psi08, psi62, psi26 and psi44: 
116  
117   a80=(-2*K80/(mu2ns*(psi100+psi82)*n))^(1/12); 
118   a08=(-2*K08/(mu2ns*(psi28+psi010)*n))^(1/12); 
119   a62=(-2*K62/(mu2ns*(psi82+psi64)*n))^(1/12); 
120   a26=(-2*K26/(mu2ns*(psi46+psi28)*n))^(1/12); 
121   a44=(-2*K44/(mu2ns*(psi64+psi46)*n))^(1/12); 
122  
123   psi80=psim('80',A,a80); 
124   psi08=psim('08',A,a08); 
125   psi62=psim('62',A,a62); 
126   psi26=psim('26',A,a26); 
127   psi44=psim('44',A,a44); 
128  
129 end; 
130  
131 if L==2 | L==3, 
132  
133   if L==2, 
134     lam80=105; 
135     lam08=105; 
136     lam62=15+90*rho^2; 
137     lam26=15+90*rho^2; 
138     lam44=9+72*rho^2+24*rho^4; 
139  
140     psi80=lam80/(64*pi*sigma1^9*sigma2*(1-rho^2)^4.5); 
141     psi08=lam08/(64*pi*sigma1*sigma2^9*(1-rho^2)^4.5); 
142     psi62=lam62/(64*pi*sigma1^7*sigma2^3*(1-rho^2)^4.5); 
143     psi26=lam26/(64*pi*sigma1^3*sigma2^7*(1-rho^2)^4.5); 
144     psi44=lam44/(64*pi*sigma1^5*sigma2^5*(1-rho^2)^4.5); 
145   end; 
146  
147 % The following are the smoothing parameters for  
148 % estimation of psi60, psi06, psi42 and psi24: 
149  
150   a60=(-2*K60/(mu2ns*(psi80+psi62)*n))^0.1; 
151   a06=(-2*K06/(mu2ns*(psi26+psi08)*n))^0.1; 
152   a42=(-2*K42/(mu2ns*(psi62+psi44)*n))^0.1; 
153   a24=(-2*K24/(mu2ns*(psi44+psi26)*n))^0.1; 
154  
155   psi60=psim('60',A,a60); 
156   psi06=psim('06',A,a06); 
157   psi42=psim('42',A,a42); 
158   psi24=psim('24',A,a24); 
159  
160 end; 
161  
162 if L==1 | L==2 | L==3, 
163  
164   if L==1, 
165     lam06=15; 
166     lam60=15; 
167     lam24=3+12*rho^2; 
168     lam42=3+12*rho^2; 
169  
170     psi60=-lam60/(32*pi*sigma1^7*sigma2*(1-rho^2)^3.5); 
171     psi06=-lam06/(32*pi*sigma1*sigma2^7*(1-rho^2)^3.5); 
172     psi42=-lam42/(32*pi*sigma1^5*sigma2^3*(1-rho^2)^3.5); 
173     psi24=-lam24/(32*pi*sigma1^3*sigma2^5*(1-rho^2)^3.5); 
174   end; 
175  
176 % The following are the smoothing parameters for  
177 % estimation of psi40, psi04 and psi22: 
178  
179   a40=(-2*K40/(mu2ns*(psi60+psi42)*n))^0.125; 
180   a04=(-2*K04/(mu2ns*(psi24+psi06)*n))^0.125; 
181   a22=(-2*K22/(mu2ns*(psi42+psi24)*n))^0.125; 
182  
183   psi40=psim('40',A,a40); 
184   psi04=psim('04',A,a04); 
185   psi22=psim('22',A,a22); 
186  
187 end; 
188  
189 if L==0, 
190  
191   lam04=3; 
192   lam40=3; 
193   lam22=3/2; 
194  
195   psi40=-lam40/(16*pi*sigma1^5*sigma2*(1-rho^2)^2.5); 
196   psi04=-lam04/(16*pi*sigma1*sigma2^5*(1-rho^2)^2.5); 
197   psi22=-lam22/(16*pi*sigma1^3*sigma2^3*(1-rho^2)^2.5); 
198  
199 end; 
200  
201  
202 h1=(psi04^0.75*R/(n*mu2^2*psi40^0.75*(sqrt(psi40*psi04)+psi22)))^(1/6); 
203 h2=h1*(psi40/psi04)^0.25; 
204  
205 h=real([h1, h2]); 
206  
207  
208 % Transform values of h back: 
209  
210 h=h*diag(std(oldA)); 
211 return 
212 function p=psim(m,A,h) 
213  
214 % PSIM        Calculate psi_m by brute force. 
215 % 
216 %             PSIM(M,A,H) 
217 % 
218 %             Christian C. Beardah 1996 
219  
220 n=length(A); 
221  
222 Xi=zeros(n,n); 
223 Xj=Xi; 
224 Yi=Xi; 
225 Yj=Xi; 
226  
227 o=ones(1,n); 
228  
229 Xj=A(:,1)*o; 
230 Xi=Xj'; 
231 Yj=A(:,2)*o; 
232 Yi=Yj'; 
233  
234 p=sum(sum(deriv2((Xi-Xj)/h,(Yi-Yj)/h,m) ))/(h*n^2); 
235 return

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