Home > wafo > kdetools > hstt.m

hstt

PURPOSE ^

Scott-Tapia-Thompson estimate of smoothing parameter.

SYNOPSIS ^

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

DESCRIPTION ^

 HSTT Scott-Tapia-Thompson estimate of smoothing parameter. 
  
  CALL: hs = hstt(data,kernel) 
  
        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 = 'epanechnikov'  - Epanechnikov kernel. (default) 
             'biweight'      - Bi-weight kernel. 
             'triweight'     - Tri-weight kernel.   
             'triangular'    - Triangular kernel. 
             'gaussian'      - Gaussian kernel 
             'rectangular'   - Rectangular kernel.  
             'laplace'       - Laplace kernel. 
             'logistic'      - Logistic kernel.   
  
  HSTT returns Scott-Tapia-Thompson (STT) estimate of smoothing 
  parameter. This is a Solve-The-Equation rule (STE). 
  Simulation studies shows that the STT estimate of HS 
  is a good choice under a variety of models. A comparison with 
  likelihood cross-validation (LCV) indicates that LCV performs slightly 
  better for short tailed densities. 
  However, STT method in contrast to LCV is insensitive to outliers. 
   
   Example:  
    x  = wnormrnd(0,1,50,1); 
    hs = hstt(x,'gauss'); 
  
  See also  hste, hbcv, hboot, hos, hldpi, hlscv, hscv, kde, kdebin

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

001 function h=hstt(A,kernel,inc,maxit,releps,abseps) 
002 %HSTT Scott-Tapia-Thompson estimate of smoothing parameter. 
003 % 
004 % CALL: hs = hstt(data,kernel) 
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 = 'epanechnikov'  - Epanechnikov kernel. (default) 
010 %            'biweight'      - Bi-weight kernel. 
011 %            'triweight'     - Tri-weight kernel.   
012 %            'triangular'    - Triangular kernel. 
013 %            'gaussian'      - Gaussian kernel 
014 %            'rectangular'   - Rectangular kernel.  
015 %            'laplace'       - Laplace kernel. 
016 %            'logistic'      - Logistic kernel.   
017 % 
018 % HSTT returns Scott-Tapia-Thompson (STT) estimate of smoothing 
019 % parameter. This is a Solve-The-Equation rule (STE). 
020 % Simulation studies shows that the STT estimate of HS 
021 % is a good choice under a variety of models. A comparison with 
022 % likelihood cross-validation (LCV) indicates that LCV performs slightly 
023 % better for short tailed densities. 
024 % However, STT method in contrast to LCV is insensitive to outliers. 
025 %  
026 %  Example:  
027 %   x  = wnormrnd(0,1,50,1); 
028 %   hs = hstt(x,'gauss'); 
029 % 
030 % See also  hste, hbcv, hboot, hos, hldpi, hlscv, hscv, kde, kdebin  
031  
032 %tested on: matlab 5.2 
033 % history 
034 % Revised pab dec2003   
035 % added inc and maxit, releps and abseps as inputs   
036 % revised pab 16.10.1999 
037 % updated string comparison to matlab 5.x 
038 % 
039 % taken from kdetools by  Christian C. Beardah 1995  
040  
041 % Reference:   
042 %  B. W. Silverman (1986)  
043 % 'Density estimation for statistics and data analysis'   
044 %  Chapman and Hall, pp 57--61  
045  
046 if nargin<2|isempty(kernel) 
047   kernel='gauss'; 
048 end 
049 if nargin<3|isempty(inc) 
050   inc = 128/2; 
051 end 
052 if nargin<4 | isempty(maxit) 
053   maxit=100; 
054 end 
055 if nargin<5 | isempty(releps) 
056   releps = 0.01; 
057 end 
058 if nargin<6 | isempty(abseps) 
059   abseps = 0.0; 
060 end 
061 [n, d] = size(A); 
062 if (n==1) & (d>1), 
063   A=A.'; 
064   n=d; 
065   d=1; 
066 end  
067  
068 [mu2,R] = kernelstats(kernel); 
069  
070 STEconstant = R /(mu2^(2)*n); 
071  
072  
073 h = hns(A,kernel); 
074  
075 % This iteration can cycle  
076 % Don't allow more than maxit iterations. 
077  
078 kopt = kdeoptset('kernel',kernel,'inc',inc); 
079  
080  
081 for dim = 1:d, 
082   count = 1; 
083   h_old = 0; 
084   h1 = h(dim); 
085    
086   while ((abs((h_old-h1))>max(releps*h1,abseps)) & (count<maxit)), 
087  
088     h_old = h1; 
089      
090     kopt.hs = h1; 
091     f = kdebin(A(:,dim),kopt); 
092  
093     delta=f.x{1}(2)-f.x{1}(1); 
094  
095     % Estimate psi4=R(f'') using simple finite differences and quadrature. 
096  
097     ix=2:(inc-1); 
098     z = ((f.f(ix+1)-2*f.f(ix)+f.f(ix-1))/delta^2).^2; 
099      
100     psi4 = delta*sum(z(:)); 
101  
102     h1 = (STEconstant/psi4)^(1/5); 
103  
104     count = count+1; 
105   end; 
106    
107   h(dim) = h1; 
108    
109   if count>= maxit 
110     disp('The obtained value did not converge.') 
111   end 
112 end

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