opt_circ_bw: Find the Optimal Bandwidth for a Circular Kernel Density Estimate

Description

This function basically wraps circular::bw.cv.ml.circular() and circular::bw.nrd.circular() of the 'circular' package, simplifying their inputs. For more control, these 'circular' functions could be used directly. The normal reference distribution ("nrd") method of finding the bandwidth parameter might give very bad results, especially for multi-modal population distributions. In these cases it can help to set kappa.est = "trigmoments".

Usage

opt_circ_bw(theta, method = c("cv", "nrd"), kappa.est = "trigmoments")

Value

A numeric value, the optimized bandwidth.

Arguments

theta: numeric vector of angles in \([-\pi, \pi)\).
method: character string describing the method, either "cv" (cross-validation), or "nrd" leading to a rule-of-thumb estimate.
kappa.est: character string describing how the spread is estimated. Either maximum likelihood "ML", or trigonometric moment "trigmoments".

Details

method="nrd" is somewhat similar to the linear case (see fit_steplength()). Instead of matching a normal distribution to the data and then calculating its optimal bandwidth, a von Mises distribution is used. To match that von Mises distribution to the data we can either find its concentration parameter kappa using maximum likelihood (kappa.est="ML") or by trigonometric moment matching (kappa.est="trigmoments"). When the data is multimodal, fitting a (unimodal) von Mises distribution using maximum likelihood will probably give bad results. Using kappa.est="trigmoments" potentially works better in those cases.

As an alternative, the bandwidth can be found by maximizing the cross-validation likelihood (method="cv"). However, with this leave-one-out cross-validation scheme, at every likelihood optimization step, \(n(n-1)\) von Mises densities need to be calculated, where \(n=\)length(theta). Therefore, this method can become quite slow with large sample sizes.

Examples

Run this code

require(circular)
require(graphics)
set.seed(123)
n <- 10  #n (number of samples) is set small for performance. Increase n to
         # a value larger than 1000 to see the effects of multimodality

angles <- rvonmisesmix(n,
  mu = c(0,pi),
  kappa = c(2,1),
  prop = c(0.5,0.5)
)
bw1 <- opt_circ_bw(theta = angles, method="nrd", kappa.est = "ML")
bw2 <- opt_circ_bw(theta = angles, method="nrd", kappa.est = "trigmoments")
bw3 <- opt_circ_bw(theta = angles, method="cv")

dens1 <- fit_angle(theta = angles, parametric = FALSE, bandwidth = bw1)
dens2 <- fit_angle(theta = angles, parametric = FALSE, bandwidth = bw2)
dens3 <- fit_angle(theta = angles, parametric = FALSE, bandwidth = bw3)
true_dens <- dvonmisesmix(
  seq(-pi,pi,0.001),
  mu = c(0,pi),
  kappa = c(2,1),
  prop = c(0.5,0.5)
)
if(interactive()){
 plot(seq(-pi, pi, 0.001), true_dens, type = "l")
 lines(as.double(dens1$x), as.double(dens1$y), col = "red")
 lines(as.double(dens2$x), as.double(dens2$y), col = "green")
 lines(as.double(dens3$x), as.double(dens3$y), col = "blue")
}