Home > wafo > kdetools > hscv.m

hscv

PURPOSE ^

Smoothed cross-validation estimate of smoothing parameter.

SYNOPSIS ^

[h,hvec,score]=hscv(A,kernel,hvec)

DESCRIPTION ^

 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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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);

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