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
High order partial derivatives of the Gaussian kernel. | |
Return 2'nd order moment of kernel pdf | |
Correlation coefficients. | |
Display message and abort function. | |
Average or mean value. | |
Standard deviation. |
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
Comments or corrections to the WAFO group