See also http://dylan-muir.com/articles/circular_kernel_estimation/
circ_ksdensityn - Compute a kernel density estimate over periodic and aperiodic domains

Usage: [vfEstimate, vfBinVol] = circ_ksdensityn(mfObservations, mfPDFSamples, <mfDomains, vfSigmas, vfWeights>)

This function calculates a kernel density estimate of an (optionally weighted) data sample, over periodic and aperiodic domains. The sample is assumed to be independent across dimensions; i.e. density estimation is performed independently for each dimension of the data.
'mfObservations' is a set of observations made over a (possibly periodic) domain. Each row corresponds to a single observation, each column corresponds to a particular dimension. By default all dimensions are periodic in [0..2*pi]; this can be modified by providing the optional argument 'mfDomains'. Each row in 'mfDomains' is [fMin fMax], one row for each dimension in 'mfObservations'. If a particular dimension should not be periodic, the corresponding row should be [nan nan]. Bounded support over a dimension is NOT implemented; each dimension is either linear and infinite or periodic.

'mfPDFSamples' defines the sample points over which to perform the kernel density estimate, over the same domains as 'mfObservations'.

Weighted estimations can be performed by providing the optional argument 'vfWeights', where each element in 'vfWeights' corresponds to the matching observation in 'mfObservations'.

The kernel density estimate will be performed using a multivariate Gaussian kernel, independent along each dimension, and wrapped along the periodic dimensions as appropriate. Kenel widths over periodic dimensions are estimated as
(4/3)^0.2 * circ_std(mfObservations(:, nDim), vfWeights) * (length(mfObservations)^-0.2)

Kernel widths over non-periodic dimensions are estimated as
(4 * std(mfObservations(:, nDim), vfWeights)^5 / 3 / length(mfObservations))^(1/5)

The optional argument 'fSigma' can be provided to set the width of the kernel.

'vfEstimate' will be a vector with a (weighted) histogram estimate of the underlying distribution, with an entry for each point in 'mfPDFSamples'. If no weighting is supplied, the estimate will be scaled to estimate a PDF over the supplied multi-dimensional domain, taking into account the estimated volume of each bin. If a weight vector is supplied, the estimate will be scaled such that the sum over the domain attempts to match the sum of weights, taking into account the multi-dimensional bin volumes.
'vfBinVol' is a vector containing volume estimates for each row in 'mfPDFSamples', under the assumption that each dimension is linearly scaled and mutually orthogonal.