circularKDE: Recent Methods for Kernel Density Estimation of Circular Data

Description

Provides recent kernel density estimation methods for circular data, including adaptive and higher-order techniques. The implementation is based on recent advances in bandwidth selection and circular smoothing. Key methods include adaptive bandwidth selection methods by Zámečník et al. (2024) doi:10.1007/s00180-023-01401-0, complete cross-validation by Hasilová et al. (2024) doi:10.59170/stattrans-2024-024, Fourier-based plug-in rules by Tenreiro (2022) doi:10.1080/10485252.2022.2057974, and higher-order kernels by Tsuruta & Sagae (2017) doi:10.1016/j.spl.2017.08.003.

Author(s)

Maintainer: Stanislav Zamecnik zamecnik@math.muni.cz

Authors:

• Ivanka Hórová horova@math.muni.cz

• Kamila Hasilová kamila.hasilova@unob.cz

• Stanislav Katina katina@math.muni.cz

See Also

Useful links:

• https://github.com/stazam/circularKDE

• Report bugs at https://github.com/stazam/circularKDE/issues

Adaptive Kernel Density Estimation for circular data

Description

This function computes an adaptive kernel density estimate for circular data, adjusting the bandwidth locally based on a global bandwidth parameter and data density (see doi:10.1007/s00180-023-01401-0).

Usage

adaptiveDensityCircular(
  x,
  bw0,
  alpha = 0.5,
  type = "n",
  z = NULL,
  from = 0,
  to = 2 * pi,
  n = 500
)

Arguments

Details

The method extends the classical idea of variable bandwidth estimators from the linear case (Breiman et al., 1977) to the circular setting by employing the von Mises kernel. Specifically, the procedure follows three steps:

1. A pilot density estimate is obtained using a fixed global bandwidth \kappa (bw0).

2. Local adaptation factors \lambda_i are computed at each data point \Theta_i according to

   \lambda_i = \left\{ g / \hat{f}(\Theta_i) \right\}^{\alpha},

   where g is a global measure of central tendency (typically the geometric mean) and \alpha \in [0,1] is a sensitivity parameter.

3. The adaptive density is evaluated using local bandwidths \lambda_i \cdot \kappa, resulting in the estimator

   \hat{f}(\theta) = \frac{1}{2n\pi} \sum_{i=1}^{n} \frac{1}{I_0(\lambda_i \kappa)} \exp\!\big(\lambda_i \kappa \cos(\theta - \Theta_i)\big).

Value

A numeric vector of length equal to the length of z (or n if z is NULL), containing the estimated density values at each evaluation point.

References

Zámečník, S., Horová, I., Katina, S., & Hasilová, K. (2023). An adaptive method for bandwidth selection in circular kernel density estimation. Computational Statistics. doi:10.1007/s00180-023-01401-0

Breiman, L., Meisel, W., & Purcell, E. (1977). Variable kernel estimates of multivariate densities. Technometrics, 19(2), 135-144. doi:10.2307/1268623

Examples

# Example with numeric data in radians
library(circular)
x <- rvonmises(100, mu = circular(0), kappa = 1)
bw0 <- bwLscvg(x = x)$minimum
dens <- adaptiveDensityCircular(x, bw0 = bw0)
plot(seq(0, 2 * pi, length.out = 500), dens, type = "l",
     main = "Adaptive Circular Density")

# Example with numerical integration over interval [0,2pi] to verify normalization of
# computed density
library(circular)
x <- rvonmises(100, mu = circular(0), kappa = 1)
bw0 <- bwLscvg(x = x)$minimum
dens <- function(z)adaptiveDensityCircular(z = z, bw0 = bw0, x=x)
integrate(dens, lower = 0, upper = 2*pi) # 1 with absolute error < 4e-07

Complete Cross-Validation for Circular Data

Description

This function calculates the optimal smoothing parameter (bandwidth) for circular data using the complete cross-validation (CCV) method (see doi:10.59170/stattrans-2024-024).

Usage

bwCcv(x, lower = 0, upper = 60, tol = 0.1)

Arguments

Details

The complete cross-validation (CCV) method is an alternative for bandwidth selection, originally proposed by Jones (1991) for linear densities. Its adaptation to the circular setting was recently studied by Hasilová et al. (2024).

The method uses functionals T_m defined as:

T_m(\kappa) = \frac{(-1)^m}{n(n-1)}\sum_{i=1}^n\sum_{j \neq i}^n K_{\kappa}^{(2m)}(\theta_{i} - \theta_{j})

where K_{\kappa}^{(2m)} is the (2m)-th derivative of K_{\kappa}.

The CCV criterion can be expressed as:

CCV(\kappa) = R(f(\kappa)) - T_0(\kappa) + \frac{1}{2}\bar{\sigma}_h^2 T_1(\kappa) + \frac{1}{24}(\eta_{2}^4(K_{\kappa}) - \eta_{4}(K_{\kappa}))T_2(\kappa)

where \eta_{j}(K_{\kappa}) denotes the j-th moment of the kernel.

For the von Mises kernel, the explicit CCV criterion becomes:

CCV(\kappa) = \frac{1}{n^2} \sum_{i=1}^n \sum_{j=1}^n (K_{\kappa} * K_{\kappa})(\theta_i - \theta_j) - T_0(\kappa) + \frac{A_1(\kappa)}{2\kappa}T_1(\kappa) + \frac{2A_1^2(\kappa) - A_2(\kappa)}{8\kappa^2}T_2(\kappa)

where A_k(\kappa) = I_k(\kappa)/I_0(\kappa) is the ratio of modified Bessel functions.

The optimal bandwidth is obtained by minimizing this criterion over the interval [lower, upper].

Value

The computed optimal smoothing parameter kappa, a numeric concentration parameter (analogous to inverse radians) that minimizes the smoothed cross-validation criterion within the interval [lower, upper] and the value of objective function at that point. Higher values indicate sharper, more concentrated kernels and less smoothing; lower values indicate broader kernels and more smoothing. If the optimization fails, a warning is issued.

References

Hasilová, K., Horová, I., Valis, D., & Zámečník, S. (2024). A comprehensive exploration of complete cross-validation for circular data. Statistics in Transition New Series, 25(3):1–12. doi:10.59170/stattrans-2024-024

See Also

bwScv, bwLscvg, bwTs

Examples

# Example with circular data
library(circular)
set.seed(123)
x <- rwrappednormal(100, mu = circular(2), rho = 0.5)
bw <- bwCcv(x)
print(round(bw$minimum, 2))

x <- rvonmises(100, mu = circular(0), kappa = 1)
bw <- bwCcv(x)
print(round(bw$minimum, 2))

Optimal Bandwidth Selection via Fourier Plug-in Method

Description

This function computes the optimal smoothing parameter (bandwidth) for circular data using the Fourier series-based direct plug-in approach based on delta sequence estimators (see doi:10.1080/10485252.2022.2057974).

Usage

bwFo(x, C1 = 0.25, C2 = 25, gamma = 0.5)

Arguments

Details

The Fourier-based plug-in estimator computes the optimal bandwidth using the formula:

\hat{h}_{FO} := (4\pi)^{-1/10} \hat{\theta}_{2,\hat{m}}^{-1/5} n^{-1/5}

where \hat{\theta}_{2,\hat{m}} is the estimator of the second-order functional \theta_2(f) based on the selected number of Fourier coefficients \hat{m}.

Under the assumption of von Mises density, this formula becomes:

\hat{h}_{VM} = (4\pi)^{-1/10} \left(\frac{3\hat{\kappa}^2 I_0(2\hat{\kappa}) - \hat{\kappa}I_1(2\hat{\kappa})}{8\pi I_0(\hat{\kappa})^2}\right)^{-1/5} n^{-1/5}

where I_0 and I_1 are the modified Bessel functions of the first kind of orders 0 and 1, and \hat{\kappa} is the estimated concentration parameter of the von Mises distribution.

Value

The computed optimal smoothing parameter kappa, a numeric concentration parameter (analogous to inverse radians) derived from the Fourier method for circular kernel density estimation. Higher values indicate sharper, more concentrated kernels and less smoothing; lower values indicate broader kernels and more smoothing.

References

Tenreiro, C. (2022). Kernel density estimation for circular data: a Fourier series-based plug-in approach for bandwidth selection. Journal of Nonparametric Statistics, 34(2):377–406. doi:10.1080/10485252.2022.2057974

See Also

bwScv, bwLscvg, bwCcv

Examples

# Example with circular data
library(circular)
set.seed(123)
x <- rvonmises(100, mu = circular(0), kappa = 2)
bw <- bwFo(x)
print(bw)

x <- rwrappednormal(100, mu = circular(1), rho = 0.7)
bw <- bwFo(x)
y <- density(x, bw=bw) 
plot(y, main="KDE with Fourier Plug-in Bandwidth")

Plug-in Method by Tsuruta and Sagae with additive Jones-Foster approach

Description

This function computes the optimal smoothing parameter (bandwidth) for circular data using the plug-in method introduced by Tsuruta and Sagae (see doi:10.1016/j.spl.2017.08.003) with the additive method from Jones and Foster (1993) to form higher-order kernel functions.

Usage

bwJf(x, verbose = FALSE)

Arguments

Details

The plug-in approach estimates the optimal bandwidth through the following steps:

1. Apply the additive method from Jones and Foster (1993) to construct a p-th order kernel function.

2. Derive expression for asymptotic mean integrated squared error (AMISE) expression.

3. Solving for the bandwidth that minimizes the AMISE. The optimal bandwidth for the additive Jones-Foster method is given by:

   \hat{\kappa}_{JF} = \left[\frac{16\sqrt{\pi}}{3} \hat{R}_{\hat{\tau}}\left(\frac{5f_{VM}^{(2)} + 2f_{VM}^{(4)}}{12}\right)n\right]^{2/9}

   where the functional \hat{R}_{\hat{\tau}} is computed as a weighted linear combination under the von Mises assumption and \hat{\tau} is the MLE estimate of the von Mises concentration parameter used as the initial value.

Value

The computed optimal smoothing parameter kappa, a numeric concentration parameter (analogous to inverse radians) derived from the circular version of the additive method for circular kernel density estimation. Higher values indicate sharper, more concentrated kernels and less smoothing; lower values indicate broader kernels and more smoothing.

References

Tsuruta, Yasuhito & Sagae, Masahiko (2017). Higher order kernel density estimation on the circle. Statistics & Probability Letters, 131:46–50. doi:10.1016/j.spl.2017.08.003

Jones, M. C. & Foster, P. J. (1993). Generalized jackknifing and higher-order kernels. Journal of Nonparametric Statistics, 3:81–94. doi:10.1080/10485259308832573

See Also

bwScv, bwLscvg, bwCcv

Examples

# Example with circular data
library(circular)
set.seed(123)
x <- rvonmises(100, mu = circular(0), kappa = 2)
bw <- bwJf(x)
print(bw)

x <- rwrappednormal(100, mu = circular(1), rho = 0.7)
bw <- bwJf(x)
print(bw)

Generalized Least Squares Cross-Validation for Circular Data

Description

This function computes the optimal smoothing parameter (bandwidth) for circular data using a generalized least squares cross-validation method (see doi:10.1007/s00180-023-01401-0).

Usage

bwLscvg(x, g = NULL, lower = 0, upper = 60, tol = 0.1)

Arguments

Details

The generalized least squares cross-validation method (LSCV_g) is an adaptation of the method originally introduced by Zhang for linear data, developed for circular data (see Zamecnik, et.al 2025) to address the finite-sample performance issues of standard LSCV.

The LSCV_g criterion is defined as:

LSCV_g(\kappa) = \frac{1}{n}R(K_{\kappa}) + \frac{1}{n(n-1)} \sum_{i=1}^n \sum_{j \neq i}^n \left(\frac{n-1}{n} (K_{\kappa}*K_{\kappa})(\theta_i-\theta_j) + \frac{2}{g(g-2)} K_{\kappa/g}(\theta_i-\theta_j) - \frac{g-1}{g-2} K_{\kappa/2}(\theta_i-\theta_j)\right)

Using the von Mises kernel, this takes the computational form:

LSCV_g(\kappa) = \frac{1}{2\pi n^2} \sum_{i=1}^n \sum_{j=1}^n \frac{I_0(\kappa \sqrt{2(1+\cos(\theta_i-\theta_j))})}{I_0^2(\kappa)} + \frac{1}{n(n-1)} \sum_{i=1}^n \sum_{j \neq i}^n \left(\frac{2}{g(g-2)} \frac{\exp(\frac{\kappa}{g}\cos(\theta_i-\theta_j))}{2\pi I_0(\kappa/g)} - \frac{g-1}{g-2} \frac{\exp(\frac{\kappa}{2}\cos(\theta_i-\theta_j))}{2\pi I_0(\kappa/2)}\right)

The optimal bandwidth is obtained by minimizing this criterion over the interval [lower, upper].

Value

The computed optimal smoothing parameter kappa, a numeric concentration parameter (analogous to inverse radians) that minimizes the smoothed cross-validation criterion within the interval [lower, upper] and the value of objective function at that point. Higher values indicate sharper, more concentrated kernels and less smoothing; lower values indicate broader kernels and more smoothing. If the optimization fails, a warning is issued.

References

Zámečník, S., Horová, I., & Hasilová, K. (2025). Generalised least square cross-validation for circular data. Communications in Statistics. doi:10.1007/s00180-023-01401-0

Zhang, J. (2015). Generalized least squares cross-validation in kernel density estimation. Statistica Neerlandica, 69(3), 315-328. doi:10.1111/stan.12061

See Also

bwScv, bwFo, bwCcv

Examples

# Example with circular data
library(circular)
set.seed(123)
x <- rwrappednormal(100, mu = circular(2), rho = 0.5)
bw <- bwLscvg(x)
print(round(bw$minimum, 2))

x <- rvonmises(100, mu = circular(0), kappa = 1)
bw <- bwLscvg(x)
print(round(bw$minimum, 2))
plot(density.circular(x, bw = bw$minimum))

Smoothed Cross-Validation for Circular Data

Description

This function computes the optimal smoothing parameter (bandwidth) for circular data using a smoothed cross-validation (SCV) method (see doi:10.1007/s00180-023-01401-0).

Usage

bwScv(x, np = 500, lower = 0, upper = 60, tol = 0.1)

Arguments

Details

The smoothed cross-validation (SCV) method is an alternative bandwidth selection approach, originally introduced by Hall & Marron (1992) for linear densities and adapted for circular data by Zámečník et al. (2023).

The SCV criterion is given by

\mathrm{SCV}(\kappa) = \frac{R(K)}{nh} + \frac{1}{n^{2}} \sum_{i=1}^{n} \sum_{j=1}^{n} \big(K_{\kappa} * K_{\kappa} * K_{\kappa} * K_{\kappa} - 2K_{\kappa} * K_{\kappa} *K_{\kappa} + K_{\kappa} * K_{\kappa}\big)(\Theta_i - \Theta_j)

where K_\kappa is the Von Mises kernel with concentration \kappa (for the formula see 3.7, 3.8 in Zámečník et al. (2023)). The optimal bandwidth minimizes the sum ISB(\kappa) + IV(\kappa) over the interval [lower, upper].

The integral expressions involved in the SCV criterion (see Sections 3.2 in Zámečník et al., 2023) are evaluated numerically using the trapezoidal rule on a uniform grid of length np.

Value

The computed optimal smoothing parameter kappa, a numeric concentration parameter (analogous to inverse radians) that minimizes the smoothed cross-validation criterion within the interval [lower, upper] and the value of objective function at that point. Higher values indicate sharper, more concentrated kernels and less smoothing; lower values indicate broader kernels and more smoothing. If the optimization fails, a warning is issued.

References

Zámečník, S., Horová, I., Katina, S., & Hasilová, K. (2023). An adaptive method for bandwidth selection in circular kernel density estimation. Computational Statistics. doi:10.1007/s00180-023-01401-0

Hall, P., & Marron, J. S. (1992). On the amount of noise inherent in bandwidth selection for a kernel density estimator. The Annals of Statistics, 20(1), 163-181.

See Also

bwTs, bwLscvg, bwCcv

Examples

# Example with circular data (Lower `nu` = more smoothing; higher = sharper peaks).
library(circular)
x <- rwrappednormal(100, mu = circular(2), rho = 0.5)
bw <- bwScv(x)
print(round(bw$minimum, 2))

x <- rvonmises(100, mu = circular(0.5), kappa = 2)
bw <- bwScv(x)
print(round(bw$minimum, 2))

Plug-in Method by Tsuruta and Sagae with multiplicative method from Terrell and Scott

Description

This function computes the optimal smoothing parameter (bandwidth) for circular data using the plug-in method, introduced by Tsuruta and Sagae (see doi:10.1016/j.spl.2017.08.003) with the multiplicative method from Terrell and Scott (1980) to form higher-order kernel functions. The method optimally balances the bias-variance tradeoff by minimizing asymptotic mean integrated squared error.

Usage

bwTs(x, verbose = FALSE)

Arguments

Details

The plug-in approach estimates the optimal bandwidth through the following steps:

1. Apply the multiplicative method from Terrell and Scott (1980) to construct a p-th order kernel function.

2. Derive expression for asymptotic mean integrated squared error (AMISE) expression.

3. Solving for the bandwidth that minimizes the AMISE. The optimal bandwidth for the multiplicative Terrell-Scott method is given by:

   \hat{\kappa}_{TS} = \left[\frac{288}{33 - 16\sqrt{2/5}} \hat{R}_{\hat{\tau}}(m_{VM}) n\right]^{2/9}

   where the computational formula is:

   m_{VM}(\theta) := [2\{f_{VM}^{(2)}\}^2/f_{VM} - 5f_{VM}^{(2)} + 2f_{VM}^{(4)}]/4

   and \hat{R}_{\hat{\tau}}(m_{VM}) is the functional computed under the von Mises assumption using the multiplicative approach. The parameter \hat{\tau} is the MLE estimate of the von Mises concentration parameter used as the initial value.

Value

The computed optimal smoothing parameter, a numeric value derived from the circular version of the multiplicative method for circular kernel density estimation.

References

Tsuruta, Yasuhito & Sagae, Masahiko (2017). Higher order kernel density estimation on the circle. Statistics & Probability Letters, 131:46–50. doi:10.1016/j.spl.2017.08.003

Terrell, George R. & Scott, David W. (1980). On improving convergence rates for nonnegative kernel density estimators. The Annals of Statistics, 8(5):1160–1163.

See Also

bwScv, bwLscvg, bwCcv

Examples

# Example with circular data
library(circular)
set.seed(123)
x <- rvonmises(100, mu = circular(0), kappa = 2)
bw <- bwTs(x)
print(bw)

x <- rwrappednormal(100, mu = circular(1), rho = 0.7)
x <- as.numeric(x)
bw <- bwTs(x)
print(bw)