HSCV Smoothed cross-validation estimate of smoothing parameter. CALL: [hs,hvec,score] = hscv(data,kernel,hvec); hs = smoothing parameter hvec = vector defining possible values of hs (default linspace(0.25*h0,h0,100), h0=hos(data,kernel)) score = score vector data = data vector kernel = 'gaussian' - Gaussian kernel the only supported Note that only the first 4 letters of the kernel name is needed. Example: data = wnormrnd(0,1,20,1) [hs hvec score] = hscv(data,'epan'); plot(hvec,score) See also hste, hbcv, hboot, hos, hldpi, hlscv, hstt, kde, kdefun
4th, 6th, 8th and 10th derivatives of the kernel function. | |
D-dimensional histogram using linear binning. | |
Oversmoothing Parameter. | |
Return 2'nd order moment of kernel pdf | |
Multivariate Kernel Function. | |
Calculates quantile levels which encloses P% of data | |
Difference and approximate derivative. | |
Discrete Fourier transform. | |
Inverse discrete Fourier transform. | |
Linearly spaced vector. | |
Standard deviation. |
001 function [h,hvec,score]=hscv(A,kernel,hvec) 002 %HSCV Smoothed cross-validation estimate of smoothing parameter. 003 % 004 % CALL: [hs,hvec,score] = hscv(data,kernel,hvec); 005 % 006 % hs = smoothing parameter 007 % hvec = vector defining possible values of hs 008 % (default linspace(0.25*h0,h0,100), h0=hos(data,kernel)) 009 % score = score vector 010 % data = data vector 011 % kernel = 'gaussian' - Gaussian kernel the only supported 012 % 013 % Note that only the first 4 letters of the kernel name is needed. 014 % 015 % Example: data = wnormrnd(0,1,20,1) 016 % [hs hvec score] = hscv(data,'epan'); 017 % plot(hvec,score) 018 % See also hste, hbcv, hboot, hos, hldpi, hlscv, hstt, kde, kdefun 019 020 021 % tested on : matlab 5.2 022 % history: 023 % revised pab Aug 2005 024 % -added support for other kernels by scaling 025 % (this is an ad hoc solution, which may be improved in future) 026 % revised pab dec2003 027 % revised pab 20.10.1999 028 % updated to matlab 5.2 029 % changed input arguments 030 % improoved gridding 031 % taken from kdetools Christian C. Beardah 1995 032 033 % Wand,M.P. and Jones, M.C. (1986) 034 % 'Kernel smoothing' 035 % Chapman and Hall, pp 75--79 036 037 % TODO % Add support for other kernels than Gaussian 038 A=A(:); 039 [n, d] = size(A); 040 if (n==1) & (d>1), 041 A=A.'; 042 n=d; 043 d=1; 044 end 045 046 if nargin<2|isempty(kernel), 047 kernel='gauss'; 048 end; 049 050 if nargin<3|isempty(hvec), 051 H=hos(A,kernel); 052 hvec=linspace(0.25*H,H,100); 053 else 054 hvec=abs(hvec); 055 end; 056 057 steps=length(hvec); 058 059 inc = 128; 060 nfft = inc*2; 061 062 xmin=min(A); % Find the minimum value of A. 063 xmax=max(A); % Find the maximum value of A. 064 xrange=xmax-xmin; % Find the range of A. 065 066 sigmaA = std(A); 067 iqr = abs(diff(qlevels2(A,[75 25]),1,1)); % interquartile range 068 k = find(iqr>0); 069 if any(k) 070 sigmaA(k) = min(sigmaA(k), iqr(k)/1.349); 071 end 072 073 % R= int(mkernel(x)^2) 074 % mu2= int(x^2*mkernel(x)) 075 [mu2,R] = kernelstats(kernel); 076 077 STEconstant = R /(mu2^(2)*n); 078 079 % xa holds the x 'axis' vector, defining a grid of x values where 080 % the k.d. function will be evaluated and plotted. 081 082 ax1=xmin-xrange/8; 083 bx1=xmax+xrange/8; 084 085 kernel2 = 'gauss'; 086 [mu2,R] = kernelstats(kernel2); 087 STEconstant2 = R /(mu2^(2)*n); 088 089 hvec = hvec*(STEconstant2/STEconstant)^(1/5); 090 for dim = 1:d 091 s = sigmaA(dim); 092 ax = ax1(dim); 093 bx = bx1(dim); 094 xa = linspace(ax,bx,inc).'; 095 xn = linspace(0,bx-ax,inc); 096 097 c = gridcount(A(:,dim),xa); 098 099 %deltax = (bx-ax)/(inc-1); 100 %xn = (0:(inc-1))*deltax; 101 %binx=floor((A(:,dim)-ax)/deltax)+1; 102 % Obtain grid counts 103 %c = full(sparse(binx,1,(xa(binx+1)-A(:,dim)),inc,1)+... 104 % sparse(binx+1,1,(A(:,dim)-xa(binx)),inc,1))/deltax; 105 106 [k40,k60,k80,k100] = deriv(0,kernel2); 107 psi8 = 105/(32*sqrt(pi)*s^9); 108 psi12 = 3465/(512*sqrt(pi)*s^13); 109 g1 = (-2*k60/(mu2*psi8*n))^(1/9); 110 %g1 = sqrt(2)*s*(2/(7*n))^(1/9) 111 g2 = (-2*k100/(mu2*psi12*n))^(1/13); 112 %g2 = sqrt(2)*s*(2/(11*n))^(1/13) 113 114 [kw4,kw6] = deriv(xn/g1,kernel2);% Obtain the kernel weights. 115 kw = [kw6,0,kw6([inc:-1:2])].';% Apply 'fftshift' to kw. 116 z = real(ifft(fft(c,nfft).*fft(kw)));% Perform the convolution. 117 118 psi6 = sum(c.*z(1:inc))/(n^2*g1^7); 119 120 % Obtain the kernel weights. 121 122 [kw4,kw6,kw8,kw10]=deriv(xn/g2,kernel2); 123 kw=[kw10,0,kw10([inc:-1:2])]';% Apply 'fftshift' to kw. 124 z=real(ifft(fft(c,nfft).*fft(kw)));% Perform the convolution. 125 126 psi10=sum(c.*z(1:inc))/(n^2*g2^11); 127 g3 = (-2*k40/(mu2*psi6*n))^(1/7); 128 % g3 = (-6/(sqrt(2*pi)*psi6*n))^(1/7) 129 g4 = (-2*k80/(mu2*psi10*n))^(1/11); 130 %g4=(-210/(sqrt(2*pi)*psi10*n))^(1/11) 131 132 % Obtain the kernel weights. 133 134 kw4=deriv(xn/g3,kernel2); 135 136 % Apply 'fftshift' to kw. 137 138 kw=[kw4,0,kw4([inc:-1:2])]'; 139 140 % Perform the convolution. 141 142 z=real(ifft(fft(c,nfft).*fft(kw))); 143 144 psi4=sum(c.*z(1:inc))/(n^2*g3^5); 145 146 % Obtain the kernel weights. 147 148 [kw4,kw6,kw8]=deriv(xn/g4,kernel2); 149 150 % Apply 'fftshift' to kw. 151 152 kw=[kw8,0,kw8([inc:-1:2])]'; 153 154 % Perform the convolution. 155 156 z=real(ifft(fft(c,nfft).*fft(kw))); 157 158 psi8=sum(c.*z(1:inc))/(n^2*g4^9); 159 160 C=(441/(64*pi))^(1/18)*(4*pi)^(-1/5)*psi4^(-2/5)*psi8^(-1/9); 161 162 M=A*ones(size(A')); 163 164 Y=(M-M'); 165 166 for i=1:steps, 167 168 g=C*n^(-23/45)*hvec(i)^(-2); 169 170 sig1=sqrt(2*hvec(i)^2+2*g^2); 171 172 sig2=sqrt(hvec(i)^2+2*g^2); 173 174 sig3=sqrt(2*g^2); 175 176 term2=sum(sum(mkernel(Y/sig1,kernel2)/sig1-2*mkernel(Y/sig2,kernel2)/sig2+mkernel(Y/sig3,kernel2)/sig3)); 177 178 score(i)=1/(n*hvec(i)*2*sqrt(pi))+term2/n^2; 179 180 end; 181 182 [L,I]=min(score); 183 184 h(dim)=hvec(I); 185 186 % Kernel other than Gaussian scale bandwidth 187 h(dim) = h(dim)*(STEconstant/STEconstant2)^(1/5); 188 end 189 190 hvec = hvec*(STEconstant/STEconstant2)^(1/5);
Comments or corrections to the WAFO group