matthiasm@8: function [bandwidth,density,xmesh]=kde(data,n,MIN,MAX)
matthiasm@8: % Reliable and extremely fast kernel density estimator for 1 dimensional data;
matthiasm@8: %        Gaussian kernel is assumed and the bandwidth is chosen automatically;
matthiasm@8: %
matthiasm@8: % INPUTS:
matthiasm@8: %     data    - a vector of data from which the density estimate is constructed;
matthiasm@8: %   MIN, MAX  - defines the interval [MIN,MAX] on which the density estimate is constructed;
matthiasm@8: %               the default values of MIN and MAX are: 
matthiasm@8: %               MIN=min(data)-Range/100 and MAX=max(data)+Range/100, where Range=max(data)-min(data);   
matthiasm@8: %          n  - the number of mesh points used in the uniform discretization of the
matthiasm@8: %               interval [MIN, MAX]; n has to be a power of two; if n is not a power of two, then
matthiasm@8: %               n is rounded up to the next power of two, i.e., n is set to n=2^ceil(log2(n));
matthiasm@8: %               the default value of n is n=2^12;
matthiasm@8: % OUTPUTS:
matthiasm@8: %   bandwidth - the optimal bandwidth (Gaussian kernel assumed);
matthiasm@8: %     density - column vector of length 'n' with the values of the density
matthiasm@8: %               estimate at the grid points;
matthiasm@8: %     xmesh   - the grid over which the density estimate is computed;
matthiasm@8: %  Reference: Botev, Z. I.,
matthiasm@8: %             "A Novel Nonparametric Density Estimator",Technical Report,The University of Queensland
matthiasm@8: %             http://espace.library.uq.edu.au/view.php?pid=UQ:12535
matthiasm@8: %
matthiasm@8: %  Example:
matthiasm@8: %    data=randn(1000,1); 
matthiasm@8: %    [bandwidth,density,xmesh]=kde(data,2^12,min(data)-1,max(data)+1);
matthiasm@8: %    plot(xmesh,density)
matthiasm@8: 
matthiasm@8: data=data(:); %make data a column vector
matthiasm@8: if nargin<2 % if n is not supplied switch to the default
matthiasm@8:     n=2^12;
matthiasm@8: end
matthiasm@8: n=2^ceil(log2(n)); % round up n to the next power of 2;
matthiasm@8: 
matthiasm@8: if nargin<4 %define the default  interval [MIN,MAX]
matthiasm@8:     minimum=min(data); maximum=max(data);
matthiasm@8:     Range=maximum-minimum;
matthiasm@8:     MIN=minimum-Range/10; MAX=maximum+Range/10;
matthiasm@8: end
matthiasm@8: % set up the grid over which the density estimate is computed;
matthiasm@8: R=MAX-MIN; dx=R/(n-1); xmesh=MIN+[0:dx:R]; N=length(data);
matthiasm@8: %bin the data uniformly using the grid define above;
matthiasm@8: initial_data=histc(data,xmesh)/N;
matthiasm@8: a=dct1d(initial_data); % discrete cosine transform of initial data
matthiasm@8: % now compute the optimal bandwidth^2 using the GCE method
matthiasm@8: t_star=gce(a,n,N);
matthiasm@8: % smooth the discrete cosine transform of initial data using t_star
matthiasm@8: a_t=a.*exp(-[0:n-1]'.^2*pi^2*t_star/2);
matthiasm@8: % now apply the inverse discrete cosine transform
matthiasm@8:   if nargout>1
matthiasm@8:     density=idct1d(a_t)/R;
matthiasm@8:   end
matthiasm@8: bandwidth=sqrt(t_star)*R;
matthiasm@8: end
matthiasm@8: function t_star=gce(a,n,N)
matthiasm@8: a=a(2:end)/2;
matthiasm@8: I=[1:n-1]'.^2;
matthiasm@8: a2=a.^2;
matthiasm@8: Var_a=zeros(n-1,1);
matthiasm@8: Var_a(1:n/2-1)=(1/2+1/2*a(2:2:n-1)-a2(1:n/2-1))/N;
matthiasm@8: 
matthiasm@8: t_star=fzero(@mise,[0,1]);
matthiasm@8: NORM=2*pi^4*sum(I.^2.*a2.*exp(-I*pi^2*t_star));
matthiasm@8: %NORM=.5*pi^4*sum([1:n-1]'.^4.*a_t(2:end).^2)/R^5;
matthiasm@8: t_star=[2*N*sqrt(pi)*NORM]^(-2/5);
matthiasm@8:     function  out=mise(t)
matthiasm@8: 
matthiasm@8:         out=sum((a2+Var_a).*(1-exp(-I*pi^2*t/2)).^2./I)+...
matthiasm@8:             sqrt(t/pi)/N*(pi^2/2)-sum(Var_a./I);
matthiasm@8:     end
matthiasm@8: end
matthiasm@8: 
matthiasm@8: function data=dct1d(data)
matthiasm@8: % computes the discrete cosine transform of the column vector data
matthiasm@8: [nrows,ncols]= size(data);
matthiasm@8: % Compute weights to multiply DFT coefficients
matthiasm@8: weight = [1;2*(exp(-i*(1:nrows-1)*pi/(2*nrows))).'];
matthiasm@8: % Re-order the elements of the columns of x
matthiasm@8: data = [ data(1:2:end,:); data(end:-2:2,:) ];
matthiasm@8: % Multiply FFT by weights:
matthiasm@8: data= real(weight.* fft(data));
matthiasm@8: end
matthiasm@8: function out = idct1d(data)
matthiasm@8: % computes the inverse discrete cosine transform
matthiasm@8: [nrows,ncols]=size(data);
matthiasm@8: % Compute weights
matthiasm@8: weights = nrows*exp(i*(0:nrows-1)*pi/(2*nrows)).';
matthiasm@8: % Compute x tilde using equation (5.93) in Jain
matthiasm@8: data = real(ifft(weights.*data));
matthiasm@8: % Re-order elements of each column according to equations (5.93) and
matthiasm@8: % (5.94) in Jain
matthiasm@8: out = zeros(nrows,1);
matthiasm@8: out(1:2:nrows) = data(1:nrows/2);
matthiasm@8: out(2:2:nrows) = data(nrows:-1:nrows/2+1);
matthiasm@8: %   Reference: 
matthiasm@8: %      A. K. Jain, "Fundamentals of Digital Image
matthiasm@8: %      Processing", pp. 150-153.
matthiasm@8: end
matthiasm@8: