Home > wafo > kdetools > hscv.m

# hscv

## PURPOSE

Smoothed cross-validation estimate of smoothing parameter.

## SYNOPSIS

[h,hvec,score]=hscv(A,kernel,hvec)

## DESCRIPTION

``` HSCV Smoothed cross-validation estimate of smoothing parameter.

CALL: [hs,hvec,score] = hscv(data,kernel,hvec);

hs     = smoothing parameter
hvec   = vector defining possible values of hs
(default linspace(0.25*h0,h0,100), h0=hos(data,kernel))
score  = score vector
data   = data vector
kernel = 'gaussian'      - Gaussian kernel the only supported

Note that only the first 4 letters of the kernel name is needed.

Example: data = wnormrnd(0,1,20,1)
[hs hvec score] = hscv(data,'epan');
plot(hvec,score)

## CROSS-REFERENCE INFORMATION

This function calls:
 deriv 4th, 6th, 8th and 10th derivatives of the kernel function. gridcount D-dimensional histogram using linear binning. hos Oversmoothing Parameter. kernelstats Return 2'nd order moment of kernel pdf mkernel Multivariate Kernel Function. 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:

## SOURCE CODE

```001 function [h,hvec,score]=hscv(A,kernel,hvec)
002 %HSCV Smoothed cross-validation estimate of smoothing parameter.
003 %
004 % CALL: [hs,hvec,score] = hscv(data,kernel,hvec);
005 %
006 %   hs     = smoothing parameter
007 %   hvec   = vector defining possible values of hs
008 %             (default linspace(0.25*h0,h0,100), h0=hos(data,kernel))
009 %   score  = score vector
010 %   data   = data vector
011 %   kernel = 'gaussian'      - Gaussian kernel the only supported
012 %
013 %  Note that only the first 4 letters of the kernel name is needed.
014 %
015 %  Example: data = wnormrnd(0,1,20,1)
016 %          [hs hvec score] = hscv(data,'epan');
017 %          plot(hvec,score)
018 % See also  hste, hbcv, hboot, hos, hldpi, hlscv, hstt, kde, kdefun
019
020
021 % tested on : matlab 5.2
022 % history:
023 % revised pab Aug 2005
024 % -added support for other kernels by scaling
025 % (this is an ad hoc solution, which may be improved in future)
026 % revised pab dec2003
027 % revised pab 20.10.1999
028 %   updated to matlab 5.2
029 % changed input arguments
030 % improoved gridding
031 % taken from kdetools     Christian C. Beardah 1995
032
033 %  Wand,M.P. and Jones, M.C. (1986)
034 % 'Kernel smoothing'
035 %  Chapman and Hall, pp 75--79
036
037 % TODO % Add support for other kernels than Gaussian
038 A=A(:);
039 [n, d] = size(A);
040 if (n==1) & (d>1),
041   A=A.';
042   n=d;
043   d=1;
044 end
045
046 if nargin<2|isempty(kernel),
047   kernel='gauss';
048 end;
049
050 if nargin<3|isempty(hvec),
051   H=hos(A,kernel);
052   hvec=linspace(0.25*H,H,100);
053 else
054   hvec=abs(hvec);
055 end;
056
057 steps=length(hvec);
058
059 inc  = 128;
060 nfft = inc*2;
061
062 xmin=min(A);    % Find the minimum value of A.
063 xmax=max(A);    % Find the maximum value of A.
064 xrange=xmax-xmin; % Find the range of A.
065
066 sigmaA = std(A);
067 iqr = abs(diff(qlevels2(A,[75 25]),1,1)); % interquartile range
068 k = find(iqr>0);
069 if any(k)
070   sigmaA(k) = min(sigmaA(k), iqr(k)/1.349);
071 end
072
073 % R= int(mkernel(x)^2)
074 % mu2= int(x^2*mkernel(x))
075 [mu2,R] = kernelstats(kernel);
076
077 STEconstant = R /(mu2^(2)*n);
078
079 % xa holds the x 'axis' vector, defining a grid of x values where
080 % the k.d. function will be evaluated and plotted.
081
082 ax1=xmin-xrange/8;
083 bx1=xmax+xrange/8;
084
085 kernel2 = 'gauss';
086 [mu2,R] = kernelstats(kernel2);
087 STEconstant2 =  R /(mu2^(2)*n);
088
089 hvec = hvec*(STEconstant2/STEconstant)^(1/5);
090 for dim = 1:d
091   s  = sigmaA(dim);
092   ax = ax1(dim);
093   bx = bx1(dim);
094   xa = linspace(ax,bx,inc).';
095   xn = linspace(0,bx-ax,inc);
096
097   c = gridcount(A(:,dim),xa);
098
099   %deltax = (bx-ax)/(inc-1);
100   %xn     = (0:(inc-1))*deltax;
101   %binx=floor((A(:,dim)-ax)/deltax)+1;
102   % Obtain grid counts
103   %c = full(sparse(binx,1,(xa(binx+1)-A(:,dim)),inc,1)+...
104   %       sparse(binx+1,1,(A(:,dim)-xa(binx)),inc,1))/deltax;
105
106   [k40,k60,k80,k100] = deriv(0,kernel2);
107   psi8  = 105/(32*sqrt(pi)*s^9);
108   psi12 = 3465/(512*sqrt(pi)*s^13);
109   g1 = (-2*k60/(mu2*psi8*n))^(1/9);
110   %g1 = sqrt(2)*s*(2/(7*n))^(1/9)
111   g2 = (-2*k100/(mu2*psi12*n))^(1/13);
112   %g2 = sqrt(2)*s*(2/(11*n))^(1/13)
113
114   [kw4,kw6] = deriv(xn/g1,kernel2);% Obtain the kernel weights.
115   kw = [kw6,0,kw6([inc:-1:2])].';% Apply 'fftshift' to kw.
116   z  = real(ifft(fft(c,nfft).*fft(kw)));% Perform the convolution.
117
118   psi6 = sum(c.*z(1:inc))/(n^2*g1^7);
119
120   % Obtain the kernel weights.
121
122   [kw4,kw6,kw8,kw10]=deriv(xn/g2,kernel2);
123   kw=[kw10,0,kw10([inc:-1:2])]';% Apply 'fftshift' to kw.
124   z=real(ifft(fft(c,nfft).*fft(kw)));% Perform the convolution.
125
126   psi10=sum(c.*z(1:inc))/(n^2*g2^11);
127   g3 = (-2*k40/(mu2*psi6*n))^(1/7);
128  % g3 = (-6/(sqrt(2*pi)*psi6*n))^(1/7)
129   g4 = (-2*k80/(mu2*psi10*n))^(1/11);
130   %g4=(-210/(sqrt(2*pi)*psi10*n))^(1/11)
131
132   % Obtain the kernel weights.
133
134   kw4=deriv(xn/g3,kernel2);
135
136   % Apply 'fftshift' to kw.
137
138   kw=[kw4,0,kw4([inc:-1:2])]';
139
140   % Perform the convolution.
141
142   z=real(ifft(fft(c,nfft).*fft(kw)));
143
144   psi4=sum(c.*z(1:inc))/(n^2*g3^5);
145
146   % Obtain the kernel weights.
147
148   [kw4,kw6,kw8]=deriv(xn/g4,kernel2);
149
150   % Apply 'fftshift' to kw.
151
152   kw=[kw8,0,kw8([inc:-1:2])]';
153
154   % Perform the convolution.
155
156   z=real(ifft(fft(c,nfft).*fft(kw)));
157
158   psi8=sum(c.*z(1:inc))/(n^2*g4^9);
159
160   C=(441/(64*pi))^(1/18)*(4*pi)^(-1/5)*psi4^(-2/5)*psi8^(-1/9);
161
162   M=A*ones(size(A'));
163
164   Y=(M-M');
165
166   for i=1:steps,
167
168     g=C*n^(-23/45)*hvec(i)^(-2);
169
170     sig1=sqrt(2*hvec(i)^2+2*g^2);
171
172     sig2=sqrt(hvec(i)^2+2*g^2);
173
174     sig3=sqrt(2*g^2);
175
176     term2=sum(sum(mkernel(Y/sig1,kernel2)/sig1-2*mkernel(Y/sig2,kernel2)/sig2+mkernel(Y/sig3,kernel2)/sig3));
177
178     score(i)=1/(n*hvec(i)*2*sqrt(pi))+term2/n^2;
179
180   end;
181
182   [L,I]=min(score);
183
184   h(dim)=hvec(I);
185
186    % Kernel other than Gaussian scale bandwidth
187   h(dim)  = h(dim)*(STEconstant/STEconstant2)^(1/5);
188 end
189
190 hvec = hvec*(STEconstant/STEconstant2)^(1/5);```

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