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

## CROSS-REFERENCE INFORMATION

This function calls:
 hns Normal Scale Estimate of Smoothing Parameter. kdebin Binned Kernel Density Estimator. kdeoptset Create or alter KDE OPTIONS structure. kernelstats Return 2'nd order moment of kernel pdf
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