Package {DirStats}

Type: Package
Title: Nonparametric Methods for Directional Data
Version: 1.0.0
Date: 2026-11-11
Description: Nonparametric kernel density estimation and bandwidth selection for directional data on the hypersphere S^q, q >= 1. Implements the kernel density estimator of Bai et al. (1988) <doi:10.1016/0047-259X(88)90113-3> in kde_dir(), the likelihood and least squares cross-validation bandwidth selectors of Hall et al. (1987) <doi:10.1093/biomet/74.4.751> in bw_dir_lcv() and bw_dir_lscv(), and the rule-of-thumb and plug-in selectors of García-Portugués (2013) <doi:10.1214/13-EJS821> in bw_dir_rot(), bw_dir_ami(), and bw_dir_emi(). Also fits mixtures of von Mises-Fisher distributions with bic_vmf_mix(), integrates on S^q, and provides fast stable evaluations of the log-Bessel function and the von Mises-Fisher normalizing constant for large orders and concentrations.
License: GPL-3
LazyData: true
Depends: R (≥ 3.6.0)
Imports: Bessel, movMF
Suggests: rotasym, testthat (≥ 3.0.0), viridisLite, withr
URL: https://github.com/egarpor/DirStats
BugReports: https://github.com/egarpor/DirStats/issues
Encoding: UTF-8
RoxygenNote: 7.3.3
Config/testthat/edition: 3
NeedsCompilation: no
Packaged: 2026-11-11 22:01:15 UTC; Eduardo
Author: Eduardo García-Portugués ORCID iD [aut, cre]
Maintainer: Eduardo García-Portugués <edgarcia@est-econ.uc3m.es>
Repository: CRAN
Date/Publication: 2026-11-12 13:00:09 UTC

DirStats – Nonparametric Methods for Directional Data

Description

Nonparametric kernel density estimation and bandwidth selection for directional data on the hypersphere \mathbb{S}^q, q \ge 1. Implements the kernel density estimator of Bai et al. (1988) doi:10.1016/0047-259X(88)90113-3 in kde_dir, the likelihood and least squares cross-validation bandwidth selectors of Hall et al. (1987) doi:10.1093/biomet/74.4.751 in bw_dir_lcv and bw_dir_lscv, and the rule-of-thumb and plug-in selectors of García-Portugués (2013) doi:10.1214/13-EJS821 in bw_dir_rot, bw_dir_ami, and bw_dir_emi. Also fits mixtures of von Mises–Fisher distributions with bic_vmf_mix, integrates on \mathbb{S}^q, and provides fast stable evaluations of the log-Bessel function and the von Mises–Fisher normalizing constant for large orders and concentrations.

Author(s)

Eduardo García-Portugués.

See Also

Useful links:

Fitting mixtures of von Mises–Fisher distributions

Description

Fitting mixtures of von Mises–Fisher distributions by the expectation–maximization algorithm, with determination of the optimal number of mixture components.

Usage

bic_vmf_mix(
  data,
  M_bound = ceiling(log(nrow(data))),
  M_neig = 3,
  crit = "BIC",
  iterative = TRUE,
  plot_it = FALSE,
  verbose = FALSE,
  kappa_max = 10000
)

Arguments

data

directional data on \mathbb{S}^q, a matrix of size c(n, q + 1).

M_bound

bound for the number of components in the mixtures. If it is not enough, the search for the mixture with minimum crit will continue from M_bound + 1 if iterative = TRUE. Defaults to ceiling(log(nrow(data))).

M_neig

number of neighbors explored around the optimal number of mixture components. Defaults to 3.

crit

information criterion to be minimized, given as the name of a function or as a function itself. It must take a fitted mixture and return a single number. Defaults to "BIC"; see the examples for a custom criterion.

iterative

keep exploring higher number of components if the optimum is attained at M_bound? Defaults to TRUE.

plot_it

display an informative plot on the optimization's grid search? Defaults to FALSE.

verbose

display fitting progress? Defaults to FALSE.

kappa_max

maximum concentration allowed in a fitted mixture. Mixtures with a component above it are discarded. Set to Inf to disable the screening. Defaults to 1e4.

Details

See Algorithm 3 in García-Portugués (2013). The expectation–maximization fit is performed with movMF.

The kappa_max screening allows discarding fits with spurious components featuring high concentration.

Value

A list with entries:

References

García-Portugués, E. (2013). Exact risk improvement of bandwidth selectors for kernel density estimation with directional data. Electronic Journal of Statistics, 7:1655–1685. doi:10.1214/13-EJS821

Hornik, K. and Grün, B. (2014). movMF: An R Package for Fitting Mixtures of von Mises–Fisher Distributions. Journal of Statistical Software, 58(10):1–31. doi:10.18637/jss.v058.i10

See Also

bw_dir_pi, vmf.

Examples

# Custom criterion: the AIC with the finite-sample correction
# 2 * k * (k + 1) / (n - k - 1), for k parameters and n observations.
AICc <- function(object) {
  ll <- logLik(object)
  k <- attr(ll, "df")
  n <- attr(ll, "nobs")
  AIC(object) + ifelse(n - k - 1 > 0, 2 * k * (k + 1) / (n - k - 1), Inf)
}

# Sample from a mixture of three von Mises--Fisher distributions
q <- 2
n <- 300
if (requireNamespace("rotasym", quietly = TRUE)) {

  set.seed(42)
  samp <- rbind(rotasym::r_vMF(n = n / 3, mu = c(rep(0, q), 1), kappa = 5),
                rotasym::r_vMF(n = n / 3, mu = c(rep(0, q), -1), kappa = 5),
                rotasym::r_vMF(n = n / 3, mu = c(1, rep(0, q)), kappa = 5))

  # Mixture fit
  print(bic_vmf_mix(data = samp, plot_it = TRUE, verbose = TRUE))
  print(bic_vmf_mix(data = samp, crit = AICc, plot_it = TRUE)$BICs)

}

Cross-validation bandwidth selectors for directional data

Description

Likelihood and least squares cross-validation bandwidth selectors for kernel density estimation with directional data.

Usage

bw_dir_lcv(
  data,
  h_grid = exp(seq(log(0.05), log(1.5), l = 100)),
  L = NULL,
  plot_it = FALSE,
  optim = TRUE,
  optim_par = 0.25,
  optim_lower = 1e-04,
  optim_upper = 10
)

bw_dir_lscv(
  data,
  h_grid = exp(seq(log(0.05), log(1.5), l = 100)),
  L = NULL,
  plot_it = FALSE,
  optim = TRUE,
  optim_par = 0.25,
  optim_lower = max(1e-04, exp(-700/(ncol(data) - 1))/sqrt(2 * pi)),
  optim_upper = 10
)

Arguments

data

directional data on \mathbb{S}^q, a matrix of size c(n, q + 1).

h_grid

vector of bandwidths for performing a grid search. Defaults to
exp(seq(log(0.05), log(1.5), l = 100)).

L

kernel function. Set internally to function(x) exp(-x) (von Mises–Fisher kernel) if NULL (default).

plot_it

display an informative plot on the optimization's grid search? Defaults to FALSE.

optim

run an optimization? Defaults to TRUE. Otherwise, a grid search on h is done. For bw_dir_lcv and bw_dir_lscv, only effective if L = NULL.

optim_par, optim_lower, optim_upper

initial value, lower bound, and upper bound of the bandwidth for optim when using the "L-BFGS-B" method. The default optim_lower of bw_dir_lscv and bw_dir_emi prevents the overflow of their losses for large q.

Details

data is not checked to have unit norm, so the user must be careful.

When L = NULL, faster vectorized expressions specific to the von Mises–Fisher kernel are employed.

The pairwise sums are accumulated in row blocks of at most 2.5e5 inner products, so the memory footprint stays flat in n for n > 5000.

bw_dir_lscv employs Monte Carlo integration for q > 2, which results in a random output. Use set.seed before to avoid it.

A warning is issued if optim does not converge, or if the selected bandwidth is at an extreme of h_grid or of c(optim_lower, optim_upper).

Value

A list with entries:

Warning

bw_dir_lscv is not suited for large dimensions, as its loss overflows for small bandwidths from q \approx 45. Use instead polykde::bw_cv_polysph with type = "LSCV" and arcsinh = TRUE.

References

Hall, P., Watson, G. S., and Cabrera, J. (1987). Kernel density estimation with spherical data. Biometrika, 74(4):751–762. doi:10.1093/biomet/74.4.751

See Also

kde_dir, bw_dir_pi.

Examples

# Sample
n <- 25
q <- 2
if (requireNamespace("rotasym", quietly = TRUE)) {

  set.seed(42)
  samp <- rotasym::r_vMF(n = n, mu = c(1, rep(0, q)), kappa = 2)

  # bw_dir_lcv
  print(bw_dir_lcv(data = samp, optim = TRUE)$h_opt)
  print(bw_dir_lcv(data = samp, optim = FALSE, plot_it = TRUE)$h_opt)
  print(bw_dir_lcv(data = samp, L = function(x) exp(-x))$h_opt)

  # bw_dir_lscv
  set.seed(42)
  print(bw_dir_lscv(data = samp, optim = TRUE)$h_opt)
  print(bw_dir_lscv(data = samp, optim = FALSE, plot_it = TRUE)$h_opt)
  print(bw_dir_lscv(data = samp, L = function(x) exp(-x))$h_opt)

}

Plug-in bandwidth selectors for directional data

Description

Plug-in bandwidth selectors for kernel density estimation with directional data, including Rule-Of-Thumb (ROT), Asymptotic MIxtures (AMI), and Exact MIxtures (EMI).

Usage

bw_dir_rot(data)

bw_dir_ami(data, fit_mix = NULL, L = NULL)

R_Psi_mixvmf(q, mu, kappa, p)

bw_dir_emi(
  data,
  fit_mix = NULL,
  optim = TRUE,
  h_grid = exp(seq(log(0.05), log(1.5), l = 100)),
  plot_it = FALSE,
  optim_par = 0.25,
  optim_lower = max(1e-04, exp(-700/(ncol(data) - 1))/sqrt(2 * pi)),
  optim_upper = 10
)

Arguments

data

directional data on \mathbb{S}^q, a matrix of size c(n, q + 1).

fit_mix

output from bic_vmf_mix. Computed internally if NULL (default).

L

kernel function. Set internally to function(x) exp(-x) (von Mises–Fisher kernel) if NULL (default).

q

dimension of \mathbb{S}^q, q \ge 1.

mu, kappa, p

mixture parameters. mu is the mean matrix of size c(length(p), q + 1), kappa is a vector of length(p) concentration parameters, and p is the vector of mixture proportions.

optim

run an optimization? Defaults to TRUE. Otherwise, a grid search on h is done. For bw_dir_lcv and bw_dir_lscv, only effective if L = NULL.

h_grid

vector of bandwidths for performing a grid search. Defaults to
exp(seq(log(0.05), log(1.5), l = 100)).

plot_it

display an informative plot on the optimization's grid search? Defaults to FALSE.

optim_par, optim_lower, optim_upper

initial value, lower bound, and upper bound of the bandwidth for optim when using the "L-BFGS-B" method. The default optim_lower of bw_dir_lscv and bw_dir_emi prevents the overflow of their losses for large q.

Details

data is not checked to have unit norm, so the user must be careful.

See Algorithms 1 (AMI) and 2 (EMI) in García-Portugués (2013). The ROT selector is implemented according to Proposition 2, but without the paper's typo in equation (6), case q = 2, where an incorrect extra \hat\kappa appears premultiplying (1 + 4 \hat\kappa^2) \sinh(2 \hat\kappa) in the denominator.

bw_dir_ami uses R_Psi_mixvmf for computing the curvature term of a mixture of von Mises–Fisher densities.

bw_dir_emi employs Monte Carlo integration for q > 2, which results in a random output. Use set.seed before to avoid it.

The selectors warn and return NA (h_opt = NA for bw_dir_emi) for degenerate concentrations.

Value

Warning

bw_dir_emi overflows for large dimensions q \approx 250, and earlier for concentrated data.

References

García-Portugués, E. (2013). Exact risk improvement of bandwidth selectors for kernel density estimation with directional data. Electronic Journal of Statistics, 7:1655–1685. doi:10.1214/13-EJS821

See Also

kde_dir, bw_dir_cv, bic_vmf_mix.

Examples

# Sample
n <- 25
q <- 2
if (requireNamespace("rotasym", quietly = TRUE)) {

  set.seed(42)
  samp <- rotasym::r_vMF(n = n, mu = c(1, rep(0, q)), kappa = 2)

  # Mixture fit
  fit_mix <- bic_vmf_mix(data = samp, plot_it = TRUE)

  # ROT
  print(bw_dir_rot(data = samp))

  # AMI
  print(bw_dir_ami(data = samp))
  print(bw_dir_ami(data = samp, fit_mix = fit_mix))
  print(bw_dir_ami(data = samp, fit_mix = fit_mix,
                   L = function(x) exp(-x)))

  # EMI
  print(bw_dir_emi(data = samp))
  print(bw_dir_emi(data = samp, fit_mix = fit_mix, optim = FALSE,
                   plot_it = TRUE))

}

Normalizing constant and moments of a directional kernel

Description

Normalizing constant c_{h, q}(L) and moments \lambda_q(L), b_q(L), and d_q(L) of a directional kernel L, which enter the kernel density estimator kde_dir and the AMISE bandwidth of bw_dir_ami.

Usage

c_h(h, q, L = NULL, log = FALSE)

lambda_L(L = NULL, q)

b_L(L = NULL, q)

d_L(L = NULL, q)

Arguments

h

vector of bandwidths.

q

dimension of \mathbb{S}^q, q \ge 1.

L

kernel function. Set internally to function(x) exp(-x) (von Mises–Fisher kernel) if NULL (default).

log

return the normalizing constant in log scale? Defaults to FALSE.

Details

The normalizing constant and the moments are given by

c_{h, q}(L)^{-1} = \omega_{q - 1} \int_{-1}^1 L\left(\frac{1 - t}{h^2}\right) (1 - t^2)^{q / 2 - 1} \,\mathrm{d}t,

\lambda_q(L) = 2^{q / 2 - 1} \omega_{q - 1} \int_0^\infty L(r) r^{q / 2 - 1} \,\mathrm{d}r,

b_q(L) = \frac{\int_0^\infty L(r) r^{q / 2} \,\mathrm{d}r} {\int_0^\infty L(r) r^{q / 2 - 1} \,\mathrm{d}r}, \quad d_q(L) = \frac{\int_0^\infty L(r)^2 r^{q / 2 - 1} \,\mathrm{d}r} {\int_0^\infty L(r) r^{q / 2 - 1} \,\mathrm{d}r},

where \omega_q is the surface area of \mathbb{S}^q (see w_q). They have closed forms for the von Mises–Fisher kernel (L = NULL), and are otherwise computed by numerical integration.

Value

References

García-Portugués, E. (2013). Exact risk improvement of bandwidth selectors for kernel density estimation with directional data. Electronic Journal of Statistics, 7:1655–1685. doi:10.1214/13-EJS821

See Also

kde_dir, bw_dir_pi.

Examples

q <- 3

# c_h
c_h(h = 0.5, q = q, L = NULL)
c_h(h = 0.5, q = q, L = function(x) exp(-x))

# b_L
b_L(L = NULL, q = q)
b_L(L = function(x) exp(-x), q = q)

# d_L
d_L(L = NULL, q = q)
d_L(L = function(x) exp(-x), q = q)

# lambda_L
lambda_L(L = NULL, q = q)
lambda_L(L = function(x) exp(-x), q = q)

Integration routines on the hypersphere

Description

Quadrature rules for integration of functions on \mathbb{S}^1, \mathbb{S}^2, and \mathbb{S}^q, q \ge 3. int_dir dispatches int_cir, int_sph, or int_hypsph, depending on q.

Usage

int_dir(f, q, sim = NULL, f_sim = NULL, ...)

int_cir(f, N = 500, na.rm = TRUE, f_vect = TRUE, ...)

int_sph(f, na.rm = TRUE, f_vect = TRUE, ...)

int_hypsph(
  f,
  q,
  M = 1e+05,
  na.rm = TRUE,
  f_vect = TRUE,
  sim = NULL,
  f_sim = NULL,
  ...
)

Arguments

f

function to be integrated on \mathbb{S}^q. Must be vectorized and accept matrix inputs of size c(nx, q + 1).

q

dimension of \mathbb{S}^q, q \ge 1.

sim

proposal for the Monte Carlo integration in int_hypsph, and in int_dir if q > 2: either a function of M returning a matrix of size c(M, q + 1), or a matrix of draws on \mathbb{S}^q. Passing a matrix reuses the same sample across several integrations. Defaults to NULL, which draws uniformly on \mathbb{S}^q.

f_sim

density of sim, used as the importance weights. Defaults to NULL, which assumes that sim is uniform.

...

further arguments passed to f. For int_dir, passed to the routine it dispatches to, and from there to f.

N

number of grid points of the trapezoidal rule. Defaults to 5e2.

na.rm

ignore possible NAs arising from the evaluation of f? Defaults to TRUE.

f_vect

can f be called in a vectorized form, with matrix input? Defaults to TRUE.

M

number of Monte Carlo replicates. Defaults to 1e5.

Details

int_cir is an extension of equation (4.1.11) in Press et al. (1997), a periodic trapezoidal rule. int_sph employs the Lebedev quadrature on \mathbb{S}^2. int_hypsph implements a Monte Carlo integration on \mathbb{S}^q. int_dir dispatches to int_cir if q = 1, to int_sph if q = 2, and to int_hypsph otherwise.

Uniform Monte Carlo misses integrands concentrated on a small region of \mathbb{S}^q, yielding estimates close to zero. For them, pass in sim a proposal that covers that region, and its density in f_sim, to perform importance-sampling Monte Carlo.

Value

A scalar approximating the integral.

References

Lebedev, V. I. and Laikov, D. N. (1999). A quadrature formula for the sphere of the 131st algebraic order of accuracy. Doklady Mathematics, 59(3):477–481.

Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P. (1997). Numerical Recipes in Fortran 77: The Art of Scientific Computing. Volume 1. Cambridge University Press, Cambridge. Second edition.

See Also

to_rad, normalize, lebedev.

Examples

# S^1, trapezoidal rule
f <- function(x) d_mixvmf(x = x, mu = c(0, 1), kappa = 2, p = 1)
int_cir(f = f)

# S^2, Lebedev rule
f <- function(x) d_mixvmf(x = x, mu = c(0, 0, 1), kappa = 2, p = 1)
int_sph(f = f)

# S^2, Monte Carlo
f <- function(x) d_mixvmf(x = x, mu = c(0, 0, 1), kappa = 2, p = 1)
int_hypsph(f = f, q = 2)

# Dispatch according to q, here to the Lebedev rule
int_dir(f = f, q = 2)

# S^2, importance-sampling Monte Carlo on a concentrated integrand, using
# a proposal that mixes the concentrated density with a uniform.
mu <- rbind(c(0, 0, 1))
kappa <- 500
f <- function(x) d_mixvmf(x = x, mu = mu, kappa = kappa, p = 1)^2
int_sph(f = f)
M_mix <- 5e3
sim <- rbind(movMF::rmovMF(n = M_mix, theta = kappa * mu, alpha = 1),
             normalize(matrix(rnorm(n = M_mix * 3), nrow = M_mix, ncol = 3)))
f_sim <- function(x) {
  0.5 * d_mixvmf(x = x, mu = mu, kappa = kappa, p = 1) + 0.5 / w_q(q = 2)
}
int_hypsph(f = f, q = 2, sim = sim, f_sim = f_sim)

Directional kernel density estimator

Description

Kernel density estimation with directional data on \mathbb{S}^q as in the estimator of Bai et al. (1988).

Usage

kde_dir(x, data, h, L = NULL)

Arguments

x

evaluation points on \mathbb{S}^q, a matrix of size c(nx, q + 1).

data

directional data on \mathbb{S}^q, a matrix of size c(n, q + 1).

h

bandwidth, a scalar.

L

kernel function. Set internally to function(x) exp(-x) (von Mises–Fisher kernel) if NULL (default).

Details

data is not checked to have unit norm, so the user must be careful.

Value

The evaluated kernel density estimator, a vector of size nx.

References

Bai, Z. D., Rao, C. R., and Zhao, L. C. (1988). Kernel estimators of density function of directional data. Journal of Multivariate Analysis, 27(1):24–39. doi:10.1016/0047-259X(88)90113-3

See Also

bw_dir_cv, bw_dir_pi, c_h, vmf.

Examples

# Sample
n <- 50
q <- 3
if (requireNamespace("rotasym", quietly = TRUE)) {

  set.seed(42)
  samp <- rotasym::r_vMF(n = n, mu = c(1, rep(0, q)), kappa = 2)

  # Evaluation points
  x <- rbind(diag(1, nrow = q + 1), diag(-1, nrow = q + 1))

  # kde_dir
  print(kde_dir(x = x, data = samp, h = 0.5, L = NULL))
  print(kde_dir(x = x, data = samp, h = 0.5, L = function(x) exp(-x)))

}

Lebedev quadrature on the sphere

Description

Nodes and weights for Lebedev quadrature on the sphere \mathbb{S}^2. The rule has 5810 points and is exact up to polynomials of order 131.

Usage

lebedev

Format

A data frame with 5810 rows and two variables:

xyz

nodes for quadrature, a matrix with three columns.

w

weights for quadrature, a vector.

Details

The approximation to the integral of f has the form

\int_{\mathbb{S}^2} f(x, y, z) \,\mathrm{d}x \,\mathrm{d}y \,\mathrm{d}z = 4 \pi \sum_{i = 1}^N w_i f(x_i, y_i, z_i)

where N = 5810. The nodes (in spherical coordinates) and weights are processed from lebedev_131.txt.

Source

https://people.sc.fsu.edu/~jburkardt/datasets/sphere_lebedev_rule/sphere_lebedev_rule.html

References

Lebedev, V. I. and Laikov, D. N. (1999). A quadrature formula for the sphere of the 131st algebraic order of accuracy. Doklady Mathematics, 59(3):477–481.

See Also

int_sph.

Examples

# Load data
data("lebedev")

# Integrate x_1 * x_2^2 (zero integral)
f_1 <- function(x) x[, 1] * x[, 2]^2
4 * pi * sum(lebedev$w * f_1(lebedev$xyz))

Stable evaluation of the scaled log-Bessel function

Description

Numerically stable evaluation of the logarithm of the exponentially scaled modified Bessel function of the first kind, \log(e^{-x} \mathcal{I}_\nu(x)). It remains finite and accurate for large orders \nu and large arguments x, where a direct call to besselI would overflow, underflow, or lose precision.

Usage

log_besselI_scaled(nu, x, spline = TRUE)

Arguments

nu

order of the Bessel function. A scalar or a vector conformable with x.

x

vector of non-negative evaluation points.

spline

evaluate the scaled log-Bessel with a fast cubic-spline interpolation? Only for the orders nu in seq(0, 50, by = 0.5). Other orders, or spline = FALSE, use besselI together with the asymptotic expansions. The spline is only effective when nu is a single value on that grid. Defaults to TRUE.

Details

The function combines four regimes so that the result is finite and accurate over the whole (\nu, x) plane:

The cutoff \nu = 50 is where besselI.nuAsym(k.max = 5) is accurate to below 10^{-10}, so it takes over from the tabulated orders with a smooth join.

The splines are built at load time from internal tables of \log(e^{-x} \mathcal{I}_\nu(x)), precomputed for the 101 orders seq(0, 50, by = 0.5) on 1400 knots equispaced in \log(x) over [10^{-5}, 10^4].

Value

The evaluation of \log(e^{-x} \mathcal{I}_\nu(x)), a vector conformable with nu and x.

References

DLMF (2026). NIST Digital Library of Mathematical Functions. Section 10.41. https://dlmf.nist.gov/10.41

See Also

log_c_vmf, besselI.

Examples

# Scaled log-Bessel across regimes
log_besselI_scaled(nu = 0.5, x = c(1e-6, 1, 1e5))

# The spline (default) and the direct evaluation agree closely
log_besselI_scaled(nu = 0.5, x = 2)
log_besselI_scaled(nu = 0.5, x = 2, spline = FALSE)

# Stable where besselI() is not
log_besselI_scaled(nu = 400, x = 20)

(Log) von Mises–Fisher normalizing constant

Description

Numerically stable evaluation of \log c_q(\kappa), the logarithm of the von Mises–Fisher normalizing constant c_q(\kappa) = \kappa^{(q - 1) / 2} / ((2\pi)^{(q + 1) / 2} \mathcal{I}_{(q - 1) / 2}(\kappa)), computed through log_besselI_scaled.

Usage

log_c_vmf(q, kappa, spline = TRUE, scaled = FALSE)

Arguments

q

dimension of \mathbb{S}^q, q\ge 0.

kappa

vector of non-negative concentration parameters.

spline

evaluate the scaled log-Bessel with a fast cubic-spline interpolation? Only for the orders nu in seq(0, 50, by = 0.5). Other orders, or spline = FALSE, use besselI together with the asymptotic expansions. The spline is only effective when nu is a single value on that grid. Defaults to TRUE.

scaled

return the exponentially scaled constant e^\kappa c_q(\kappa)? Defaults to FALSE.

Details

The value at the uniform limit \kappa = 0 is the reciprocal surface area of \mathbb{S}^q, obtained by continuity.

Value

The evaluation of \log c_q(\kappa), a vector conformable with q and kappa.

See Also

log_besselI_scaled, vmf.

Examples

# Normalizing constant, uniform case included
log_c_vmf(q = 2, kappa = c(0, 1, 10))

# Agrees with rotasym::c_vMF() for moderate concentrations
if (requireNamespace("rotasym", quietly = TRUE)) {

  rotasym::c_vMF(p = 3, kappa = c(0, 1, 10), log = TRUE)

}

Normalization onto the hypersphere

Description

Euclidean norms of the rows of x, and normalization of data in \mathbb{R}^{q + 1} to \mathbb{S}^q.

Usage

norm2(x)

normalize(x)

Arguments

x

a matrix of size c(n, q + 1) or a vector of length q + 1.

Value

See Also

int_dir.

Examples

# Normalization
x <- 1:3
norm2(x)
normalize(x)
x <- rbind(1:3, 3:1)
norm2(x)
normalize(x)

Polar and spherical coordinates

Description

Transformations between \mathbb{S}^1 and [0, 2\pi), and between \mathbb{S}^2 and [0, 2\pi) \times [0, \pi].

Usage

to_cir(th)

to_rad(x)

to_sph(th, ph)

Arguments

th

vector of angles in [0, 2\pi).

x

points on \mathbb{S}^1 or \mathbb{S}^2, a matrix of size c(n, 2) or c(n, 3), or a vector.

ph

vector of angles in [0, \pi].

Value

See Also

int_dir.

Examples

# Circular transformations
th <- 1
x <- c(0, 1)
to_rad(to_cir(th))
to_rad(to_cir(c(th, th + 1)))
to_cir(to_rad(x))
to_cir(to_rad(rbind(x, -x)))

# Spherical transformations
th <- 2
ph <- 1
x <- c(0, 1, 0)
to_rad(to_sph(th, ph))
to_rad(to_sph(c(th, th + 1),
              c(ph, ph + 1)))
to_sph(to_rad(x)[, 1], to_rad(x)[, 2])
to_sph(to_rad(rbind(x, -x))[, 1], to_rad(rbind(x, -x))[, 2])

Von Mises–Fisher distribution utilities

Description

Maximum likelihood estimation for the von Mises–Fisher distribution and evaluation of density mixtures.

Usage

kappa_ml(data, min_kappa = 1e-04, max_kappa = 10000, R = NULL, ...)

mu_ml(data)

d_mixvmf(x, mu, kappa, p, norm = FALSE)

Arguments

data

directional data on \mathbb{S}^q, a matrix of size c(n, q + 1).

min_kappa, max_kappa

initial search interval (in kappa) for the maximum likelihood estimate. The interval is automatically extended if the root lies outside. Defaults to 1e-4 and 1e4.

R

resultant length, a scalar. If NULL (default), it is computed from data.

...

further parameters passed to uniroot.

x

evaluation points on \mathbb{S}^q, a matrix of size c(nx, q + 1).

mu, kappa, p

mixture parameters. mu is the mean matrix of size c(length(p), q + 1), kappa is a vector of length(p) concentration parameters, and p is the vector of mixture proportions.

norm

enforce normalization of x internally? Defaults to FALSE.

Value

See Also

kde_dir, bic_vmf_mix.

Examples

# Sample
n <- 50
q <- 2
if (requireNamespace("rotasym", quietly = TRUE)) {

  set.seed(42)
  samp <- rotasym::r_vMF(n = n, mu = c(1, rep(0, q)), kappa = 2)

  # Estimates
  print(mu_ml(data = samp))
  print(kappa_ml(data = samp))

}

# Mixture
x <- to_cir(seq(0, 2 * pi, l = 200))
dens <- d_mixvmf(x = x, mu = rbind(c(-1, 0), c(0, 1), c(1, 0)),
                 kappa = 1:3, p = c(0.5, 0.2, 0.3))
plot(to_rad(x), dens, type = "l")

Surface area of the hypersphere

Description

Surface area of the hypersphere \mathbb{S}^q := \{x \in \mathbb{R}^{q + 1} : ||x|| = 1\}, q \ge 0, given by \omega_q := 2\pi^{(q + 1) / 2} / \Gamma((q + 1) / 2).

Usage

w_q(q, log = FALSE)

Arguments

q

dimension of \mathbb{S}^q, q \ge 0. A vector.

log

return the logarithm of the area? Defaults to FALSE.

Details

\mathbb{S}^0 = \{-1, 1\} and its "surface area" is 2. The function is vectorized on q.

Value

The surface area of \mathbb{S}^q, a vector of size length(q).

See Also

int_hypsph, log_c_vmf.

Examples

# Areas of S^0, S^1, and S^2
w_q(q = 0:2)

# Area as a function of q
q <- 0:19
plot(q, w_q(q = q), type = "o", pch = 16, xlab = "q", ylab = "Area",
     main = expression("Surface area of " * S^q))