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'); See also hbcv, hboot, hos, hldpi, hlscv, hscv, hstt, kde, kdefun
4th, 6th, 8th and 10th derivatives of the kernel function. | |
D-dimensional histogram using linear binning. | |
Normal Scale Estimate of Smoothing Parameter. | |
Return 2'nd order moment of kernel pdf | |
Calculates quantile levels which encloses P% of data | |
Difference and approximate derivative. | |
Discrete Fourier transform. | |
Inverse discrete Fourier transform. | |
Linearly spaced vector. | |
Standard deviation. |
Biased Cross-Validation smoothing parameter for 2D data. |
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 032 % added todo comments 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
Comments or corrections to the WAFO group