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

## CROSS-REFERENCE INFORMATION

This function calls:
 deriv 4th, 6th, 8th and 10th derivatives of the kernel function. gridcount D-dimensional histogram using linear binning. hns Normal Scale Estimate of Smoothing Parameter. kernelstats Return 2'nd order moment of kernel pdf qlevels2 Calculates quantile levels which encloses P% of data diff Difference and approximate derivative. fft Discrete Fourier transform. ifft Inverse discrete Fourier transform. linspace Linearly spaced vector. std Standard deviation.
This function is called by:
 hbcv2 Biased Cross-Validation smoothing parameter for 2D data.

## 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
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