Home > wafo > kdetools > hste.m

hste

PURPOSE ^

2-Stage Solve the Equation estimate of smoothing parameter.

SYNOPSIS ^

h=hste(A,kernel,h,inc,maxit,releps,abseps)

DESCRIPTION ^

 HSTE 2-Stage Solve the Equation estimate of smoothing parameter. 
  
  CALL:  hs = hste(data,kernel,h0) 
   
        hs = one dimensional value for smoothing parameter 
             given the data and kernel.  size 1 x D 
    data   = data matrix, size N x D (D = # dimensions ) 
    kernel = 'gaussian'  - Gaussian kernel (default) 
              ( currently the only supported kernel) 
        h0 = initial starting guess for hs (default h0=hns(A,kernel)) 
  
   Example:  
    x  = wnormrnd(0,1,50,1); 
    hs = hste(x,'gauss'); 
  
  See also  hbcv, hboot, hos, hldpi, hlscv, hscv, hstt, kde, kdefun

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function h=hste(A,kernel,h,inc,maxit,releps,abseps) 
002 %HSTE 2-Stage Solve the Equation estimate of smoothing parameter. 
003 % 
004 % CALL:  hs = hste(data,kernel,h0) 
005 %  
006 %       hs = one dimensional value for smoothing parameter 
007 %            given the data and kernel.  size 1 x D 
008 %   data   = data matrix, size N x D (D = # dimensions ) 
009 %   kernel = 'gaussian'  - Gaussian kernel (default) 
010 %             ( currently the only supported kernel) 
011 %       h0 = initial starting guess for hs (default h0=hns(A,kernel)) 
012 % 
013 %  Example:  
014 %   x  = wnormrnd(0,1,50,1); 
015 %   hs = hste(x,'gauss'); 
016 % 
017 % See also  hbcv, hboot, hos, hldpi, hlscv, hscv, hstt, kde, kdefun 
018  
019 % Reference:   
020 %  B. W. Silverman (1986)  
021 % 'Density estimation for statistics and data analysis'   
022 %  Chapman and Hall, pp 57--61 
023 % 
024 %  Wand,M.P. and Jones, M.C. (1986)  
025 % 'Kernel smoothing' 
026 %  Chapman and Hall, pp 74--75 
027    
028 % tested on: matlab 5.2 
029 % revised pab aug 2005 
030 % - All kernels supported 
031 % Revised pab dec 2003 
032 % added todo comments   
033 % added inc, maxit,abseps,releps as inputs   
034 % revised pab 16.10.1999 
035 % added h0 as input 
036 %  the gridding is made much faster 
037 % taken from kdetools   Christian C. Beardah 1995  
038  
039 % TODO % NB: this routine can be made faster: 
040 % TODO % replace the iteration in the end with a Newton Raphson scheme 
041  
042    
043  
044  
045  
046 if nargin<2|isempty(kernel) 
047  kernel='gauss'; 
048 end 
049 if nargin<3|isempty(h) 
050   h=hns(A,kernel); 
051 end 
052 if nargin<4 | isempty(inc) 
053   inc=128; 
054 end 
055 if nargin<5 | isempty(maxit) 
056   maxit = 100; 
057 end 
058 if nargin<6 | isempty(releps) 
059   releps = 0.01; 
060 end 
061 if nargin<7 | isempty(abseps) 
062   abseps = 0.0; 
063 end 
064 [n, d] = size(A); 
065 if (n==1) & (d>1), 
066   A=A.'; 
067   n=d; 
068   d=1; 
069 end 
070  
071  
072 % R   = int(mkernel(x)^2) 
073 % mu2 = int(x^2*mkernel(x)) 
074 [mu2,R] = kernelstats(kernel); 
075  
076 STEconstant = R /(mu2^(2)*n); 
077  
078  
079 nfft = inc*2; 
080  
081  
082 xmin   = min(A);    % Find the minimum value of A. 
083 xmax   = max(A);    % Find the maximum value of A. 
084 xrange = xmax-xmin; % Find the range of A. 
085  
086 sigmaA = std(A); 
087 iqr = abs(diff(qlevels2(A,[75 25]),1,1)); % interquartile range 
088 k = find(iqr>0); 
089 if any(k) 
090   sigmaA(k) = min(sigmaA(k), iqr(k)/1.349); 
091 end 
092  
093 % xa holds the x 'axis' vector, defining a grid of x values where  
094 % the k.d. function will be evaluated. 
095  
096 ax1 = xmin-xrange/8; 
097 bx1 = xmax+xrange/8; 
098  
099 kernel2 = 'gaus'; % kernel 
100 [mu2,R] = kernelstats(kernel2); 
101 STEconstant2 = R /(mu2^(2)*n); 
102 for dim = 1:d 
103   s = sigmaA(dim); 
104   ax = ax1(dim); 
105   bx = bx1(dim); 
106    
107   xa = linspace(ax,bx,inc).';  
108   xn = linspace(0,bx-ax,inc); 
109    
110   c = gridcount(A(:,dim),xa); 
111    
112   %deltax = (bx-ax)/(inc-1); 
113   %xn     = (0:(inc-1))*deltax; 
114   %binx   = floor((A(:,dim)-ax)/deltax)+1; 
115   % Obtain grid counts 
116   %c = full(sparse(binx,1,(xa(binx+1)-A(:,dim)),inc,1)+... 
117     %   sparse(binx+1,1,(A(:,dim)-xa(binx)),inc,1))/deltax; 
118  
119  
120   % Step 1 
121   psi6NS = -15/(16*sqrt(pi)*s^7); 
122   psi8NS = 105/(32*sqrt(pi)*s^9); 
123  
124   % Step 2 
125   [k40,k60] = deriv(0,kernel2); 
126   g1 = (-2*k40/(mu2*psi6NS*n))^(1/7); 
127   g2 = (-2*k60/(mu2*psi8NS*n))^(1/9); 
128  
129   % Estimate psi6 given g2. 
130   [kw4,kw6] = deriv(xn/g2,kernel2); % kernel weights. 
131   kw   = [kw6,0,kw6([inc:-1:2])].';             % Apply 'fftshift' to kw. 
132   z    = real(ifft(fft(c,nfft).*fft(kw)));     % convolution. 
133   psi6 = sum(c.*z(1:inc))/(n*(n-1)*g2^7); 
134  
135   % Estimate psi4 given g1. 
136   kw4  = deriv(xn/g1,kernel2); % kernel weights. 
137   kw   = [kw4,0,kw4([inc:-1:2])]';            % Apply 'fftshift' to kw. 
138   z    = real(ifft(fft(c,nfft).*fft(kw)));    % convolution. 
139   psi4 = sum(c.*z(1:inc))/(n*(n-1)*g1^5); 
140  
141  
142  
143   % 
144   h1    = h(dim); 
145   h_old = 0; 
146   count = 0; 
147    
148   while ((abs(h_old-h1)>max(releps*h1,abseps)) & (count < maxit)), 
149     count = count+1; 
150     % save old value 
151     h_old = h1; 
152    
153     % Step 3 
154     gamma=((2*k40*mu2*psi4*h1^5)/(-psi6*R))^(1/7); 
155  
156     % Now estimate psi4 given gamma. 
157     kw4 = deriv(xn/gamma,kernel2); %kernel weights.  
158     kw  = [kw4,0,kw4([inc:-1:2])]'; % Apply 'fftshift' to kw. 
159     z   = real(ifft(fft(c,nfft).*fft(kw))); % convolution. 
160  
161     psi4Gamma  = sum(c.*z(1:inc))/(n*(n-1)*gamma^5); 
162    
163     % Step 4 
164     h1 = (STEconstant2/psi4Gamma)^(1/5); 
165   end;   
166   % Kernel other than Gaussian scale bandwidth 
167   h1  = h1*(STEconstant/STEconstant2)^(1/5); 
168    
169  
170   if count>= maxit 
171     disp('The obtained value did not converge.') 
172   end 
173   h(dim) = h1; 
174 end % for dim loop

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