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
Normal Scale Estimate of Smoothing Parameter. | |
Binned Kernel Density Estimator. | |
Create or alter KDE OPTIONS structure. | |
Return 2'nd order moment of kernel pdf |
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
Comments or corrections to the WAFO group