Code covered by the BSD License  

Highlights from
Kernel density estimation for circular functions

image thumbnail

Kernel density estimation for circular functions

by

Dylan Muir (view profile)

 

24 Oct 2013 (Updated )

Performs kernel density estimates over arbitrary periodic domains.

circ_ksdensity(vfObservations, vfPDFSamples, vfDomain, fSigma, vfWeights)
function [vfEstimate] = circ_ksdensity(vfObservations, vfPDFSamples, vfDomain, fSigma, vfWeights)

% circ_ksdensity FUNCTION - Compute a kernel density estimate over a periodic domain
%
% Usage: [vfEstimate] = circ_ksdensity(vfObservations, vfPDFSamples, 
%                                      <vfDomain, fSigma, vfWeights>)
%
% This function calculates a kernel density estimate of an (optionally
% weighted) data sample, over a periodic domain.
%
% 'vfObservations' is a set of observations made over a periodic domain,
% optionally defined by 'vfDomain': [fMin fMax]. The default domain is
% [0..2*pi]. 'vfPDFSamples' defines the sample points over which to perform
% the kernel density estimate, over the same domain as 'vfObservations'.
%
% Weighted estimations can be performed by providing the optional argument
% 'vfWeights', where each element in 'vfWeights' corresponds to the
% matching element in 'vfObservations'.
%
% The kernel density estimate will be performed using a wrapped Gaussian
% kernel, with a width estimated as
%    (4/3)^0.2 * circ_std(vfObservations, vfWeights) * (length(vfObservations^-0.2)
%
% The optional argument 'fSigma' can be provided to set the width of the
% kernel.
%
% 'vfEstimate' will be a vector with a (weighted) estimate of the
% underlying distribution, with an entry for each element of
% 'vfPDFSamples'. If no weighting is supplied, the estimate will be scaled
% such that it forms a PDF estimate over the supplied sample domain, taking
% into account sample bin widths. If a weight vector is supplied then the
% estimate will be scaled such that the sum over the domain attempts to
% match the sum of weights, taking into account sample bin widths.

% Author: Dylan Muir <dylan.muir@unibas.ch>
% Created: 23rd October, 2013

% -- Defaults

DEF_vfDomain = [0 2*pi];


% -- Check arguments

if (nargin < 2)
   help circ_ksdensity;
   error('circ_ksdensity:Usage', ...
         '*** circ_ksdensity: Incorrect usage');
end

if (~exist('vfDomain', 'var') || isempty(vfDomain))
   vfDomain = DEF_vfDomain;
end

% - Do we need to estimate fSigma?
if (~exist('fSigma', 'var'))
   % - Sigma will be estimated
   fSigma = [];
end

vfObservations = vfObservations(:);
vnPSFSamplesSize = size(vfPDFSamples);
vfPDFSamples = vfPDFSamples(:);

% - If weights are not provided, weight each observation equally
if (~exist('vfWeights', 'var'))
   vfWeights = ones(size(vfObservations)) ./ numel(vfObservations);

% - Check the number of observations matches the number of weights
elseif (numel(vfObservations) ~= numel(vfWeights))
   error('circ_ksdensity:Usage', ...
         '*** circ_ksdensity: The number of observations must be equal to the number of weights.');
end


% -- Map everything to [0..2pi] and wrap over domain

vfObservations = (vfObservations - vfDomain(1)) ./ diff(vfDomain) .* 2*pi;
vfObservations = mod(vfObservations, 2*pi);

vfPDFSamples = (vfPDFSamples - vfDomain(1)) ./ diff(vfDomain) .* 2*pi;
vfPDFSamples = mod(vfPDFSamples, 2*pi);


% - Estimate sigma, if necessary
if (isempty(fSigma))
   fSigma = (4/3)^0.2 * circ_std(vfObservations, vfWeights) * (numel(vfObservations)^-0.2);
end


% -- Pad observations above and below domain

vfObservations = [vfObservations;
                  vfObservations - 2*pi;
                  vfObservations + 2*pi];
vfWeights = repmat(vfWeights, 3, 1) ./ 3;


% -- Perform kernel density estimate

vfEstimate = ksdensity(vfObservations, vfPDFSamples, 'weights', vfWeights', 'width', fSigma);

% - Reshape return to match shape of 'vfPDFSamples'
vfEstimate = reshape(vfEstimate, vnPSFSamplesSize);

% - Correct scaling of histogram estimate
vfEstimate = vfEstimate .* 3 .* sum(vfWeights);

% --- END of circ_ksdensity FUNCTION ---


% circ_std FUNCTION Estimate the weighted circular standard deviation of a dataset
function [s s0] = circ_std(alpha, w)

% compute mean resultant vector length
r = circ_r(alpha,w);

s = sqrt(2*(1-r));      % 26.20
s0 = sqrt(-2*log(r));    % 26.21


% circ_r FUNCTION Compute the weighted resultant of a dataset
function r = circ_r(alpha, w)

% compute weighted sum of cos and sin of angles
r = sum(w.*exp(1i*alpha),1);

% obtain length
r = abs(r)./sum(w,1);


% --- END of circ_ksdensity.m ---


Contact us