Package 'polykde'

Description

Kernel density estimation on the polysphere, hypersphere, and circle. Includes functions for density estimation, regression estimation, ridge estimation, bandwidth selection, kernels, samplers, and homogeneity tests. Companion package to García-Portugués and Meilán-Vila (2024) <doi:10.48550/arXiv.2411.04166> and García-Portugués and Meilán-Vila (2023) <doi:10.1007/978-3-031-32729-2_4>.

Author(s)

Eduardo García-Portugués.

References

García-Portugués, E. and Meilán-Vila, A. (2024). Kernel density estimation with polyspherical data and its applications. arXiv:2411.04166. doi:10.48550/arXiv.2411.04166.

García-Portugués, E. and Meilán-Vila, A. (2023). Hippocampus shape analysis via skeletal models and kernel smoothing. In Larriba, Y. (Ed.), Statistical Methods at the Forefront of Biomedical Advances, pp. 63–82. Springer, Cham. doi:10.1007/978-3-031-32729-2_4.

See Also

Useful links:

• https://github.com/egarpor/polykde

• Report bugs at https://github.com/egarpor/polykde/issues

Description

Obtain the angular coordinates of points on a polysphere $\mathcal{S}^{d_1}\times\cdots\times\mathcal{S}^{d_r}$, and vice versa.

Usage

angles_to_polysph(theta, d)

polysph_to_angles(x, d)

Arguments

Value

• angles_to_polysph: the matrix x.

• polysph_to_angles: the matrix theta.

Examples

# Check changes of coordinates
polysph_to_angles(angles_to_polysph(rep(pi / 2, 3), d = 2:1), d = 2:1)
angles_to_polysph(polysph_to_angles(x = c(0, 0, 1, 0, 1), d = 2:1), d = 2:1)

Description

Transforms the angles $(\theta_1,\ldots,\theta_d)$ in $[0,\pi)^{d-1}\times[-\pi,\pi)$ into the Cartesian coordinates

$(\cos(x_1),\sin(x_1)\cos(x_2),\ldots, \sin(x_1)\cdots\sin(x_{d-1})\cos(x_d), \sin(x_1)\cdots\sin(x_{d-1})\sin(x_d))$

of the hypersphere $\mathcal{S}^{d}$, and vice versa.

Usage

angles_to_sph(theta)

sph_to_angles(x)

Arguments

Value

• angles_to_sph: the matrix x.

• sph_to_angles: the matrix theta.

Examples

# Check changes of coordinates
sph_to_angles(angles_to_sph(c(pi / 2, 0, pi)))
sph_to_angles(angles_to_sph(rbind(c(pi / 2, 0, pi), c(pi, pi / 2, 0))))
angles_to_sph(sph_to_angles(c(0, sqrt(0.5), sqrt(0.1), sqrt(0.4))))
angles_to_sph(sph_to_angles(rbind(c(0, sqrt(0.5), sqrt(0.1), sqrt(0.4)),
                                  c(0, sqrt(0.5), sqrt(0.5), 0),
                                  c(0, 1, 0, 0),
                                  c(0, 0, 0, -1),
                                  c(0, 0, 1, 0))))

Description

Transforms the angles $(\theta_1,\ldots,\theta_d)$ in $[-\pi,\pi)^d$ into the Cartesian coordinates

$(\cos(x_1), \sin(x_1),\ldots,\cos(x_d), \sin(x_d))$

of the torus $(\mathcal{S}^1)^d$, and vice versa.

Usage

angles_to_torus(theta)

torus_to_angles(x)

Arguments

Value

• angles_to_torus: the matrix x.

• torus_to_angles: the matrix theta.

Examples

# Check changes of coordinates
torus_to_angles(angles_to_torus(c(0, pi / 3, pi / 2)))
torus_to_angles(angles_to_torus(rbind(c(0, pi / 3, pi / 2), c(0, 1, -2))))
angles_to_torus(torus_to_angles(c(0, 1, 1, 0)))
angles_to_torus(torus_to_angles(rbind(c(0, 1, 1, 0), c(0, 1, 0, 1))))

Description

Computes least squares cross-validation bandwidths for kernel regression estimation with polyspherical response and scalar predictor. It computes both the bandwidth that minimizes the cross-validation loss and its "one standard error" variant.

Usage

bw_cv_kre_polysph(X, Y, d, p = 0, h_grid = bw.nrd(X) * 10^seq(-2, 2, l =
  100), plot_cv = TRUE, fast = TRUE)

Arguments

Details

A similar output to glmnet's cv.glmnet is returned.

Value

A list with the following fields:

Examples

n <- 50
X <- seq(0, 1, l = n)
Y <- r_path_s2r(n = n, r = 1, sigma = 0.1, spiral = TRUE)[, , 1]
bw_cv_kre_polysph(X = X, Y = Y, d = 2, p = 0)
bw_cv_kre_polysph(X = X, Y = Y, d = 2, p = 1, fast = FALSE)

Description

Likelihood Cross-Validation (LCV) and Least Squares Cross-Validation (LSCV) bandwidth selection for the polyspherical kernel density estimator.

Usage

bw_cv_polysph(X, d, kernel = 1, kernel_type = 1, k = 10,
  type = c("LCV", "LSCV")[1], M = 10000, bw0 = NULL, na.rm = FALSE,
  ncores = 1, h_min = 0, upscale = FALSE, deriv = 0, ...)

Arguments

Value

A list as optim or optimParallel output. In particular, the optimal bandwidth is stored in par.

Examples

n <- 50
d <- 1:2
kappa <- rep(10, 2)
X <- r_vmf_polysph(n = n, d = d, mu = r_unif_polysph(n = 1, d = d),
                   kappa = kappa)
bw_cv_polysph(X = X, d = d, type = "LCV")$par
bw_cv_polysph(X = X, d = d, type = "LSCV")$par

Description

This function computes the minimum bandwidth allowed in likelihood cross-validation with Epanechnikov kernels, for a given dataset and dimension.

Usage

bw_lcv_min_epa(X, d, kernel_type = c("prod", "sph")[1])

Arguments

Value

The minimum bandwidth allowed.

Examples

n <- 5
d <- 1:3
X <- r_unif_polysph(n = n, d = d)
h_min <- rep(bw_lcv_min_epa(X = X, d = d), length(d))
log_cv_kde_polysph(X = X, d = d, h = h_min - 1e-4, kernel = 2) # Problem
log_cv_kde_polysph(X = X, d = d, h = h_min + 1e-4, kernel = 2) # OK

Description

Computes marginal (sphere by sphere) rule-of-thumb bandwidths for the polyspherical kernel density estimator using a von Mises–Fisher distribution as reference.

Usage

bw_mrot_polysph(X, d, kernel = 1, k = 10, upscale = FALSE, deriv = 0,
  kappa = NULL, ...)

Arguments

Value

A vector of size r with the marginal optimal bandwidths.

Examples

n <- 100
d <- 1:2
kappa <- rep(10, 2)
X <- r_vmf_polysph(n = n, d = d, mu = r_unif_polysph(n = 1, d = d),
                   kappa = kappa)
bw_rot_polysph(X = X, d = d)$bw
bw_mrot_polysph(X = X, d = d)

Description

Computes the rule-of-thumb bandwidth for the polyspherical kernel density estimator using a product of von Mises–Fisher distributions as reference in the Asymptotic Mean Integrated Squared Error (AMISE).

Usage

bw_rot_polysph(X, d, kernel = 1, kernel_type = c("prod", "sph")[1],
  bw0 = NULL, upscale = FALSE, deriv = 0, k = 10, kappa = NULL, ...)

Arguments

Details

The selector assumes that the density curvature matrix $\boldsymbol{R}$ of the unknown density is approximable by that of a product of von Mises–Fisher densities, $\boldsymbol{R}(\boldsymbol{\kappa})$. The estimation of the concentration parameters $\boldsymbol{\kappa}$ is done by maximum likelihood.

Value

A list with entries bw (optimal bandwidth) and opt, the latter containing the output of nlm.

Examples

n <- 100
d <- 1:2
kappa <- rep(10, 2)
X <- r_vmf_polysph(n = n, d = d, mu = r_unif_polysph(n = 1, d = d),
                   kappa = kappa)
bw_rot_polysph(X = X, d = d)$bw

Description

Remove points from the ridge that are spurious. The cleaning is done by removing end points in the Euler algorithm that did not converge, do not have a negative second eigenvalue, or are in low-density regions.

Usage

clean_euler_ridge(e, X, p_out = NULL)

Arguments

Value

A list with the same structure as that returned by euler_ridge, but with the spurious points. The removed points are informed in the removed field.

Examples

## Test on S^2 with some spurious end points

# Sample
r <- 1
d <- 2
n <- 50
ind_dj <- comp_ind_dj(d = d)
set.seed(987202226)
X <- r_path_s2r(n = n, r = r, spiral = FALSE, Theta = cbind(c(1, 0, 0)),
                sigma = 0.2)[, , 1]
col_X <- rep(gray(0), n)
col_X_alp <- rep(gray(0, alpha = 0.25), n)

# Euler
h_rid <- 0.5
h_eu <- h_rid^2
N <- 30
eps <- 1e-6
X0 <- r_unif_polysph(n = n, d = d)
Y <- euler_ridge(x = X0, X = X, d = d, h = h_rid, h_euler = h_eu,
                 N = N, eps = eps, keep_paths = TRUE)
Y_removed <- clean_euler_ridge(e = Y, X = X)$removed
col_X[Y_removed] <- 2
col_X_alp[Y_removed] <- 2

# Visualization
i <- N # Between 1 and N
sc3 <- scatterplot3d::scatterplot3d(Y$paths[, , 1], color = col_X_alp,
                                    pch = 19, xlim = c(-1, 1),
                                    ylim = c(-1, 1), zlim = c(-1, 1),
                                    xlab = "x", ylab = "y", zlab = "z")
sc3$points3d(rbind(Y$paths[, , i]), col = col_X, pch = 16, cex = 0.75)
for (k in seq_len(nrow(Y$paths))) {

  sc3$points3d(t(Y$paths[k, , ]), col = col_X_alp[k], type = "l")

}

Description

Given Cartesian coordinates of polyspherical data, computes the 0-based indexes at which the Cartesian coordinates for each hypersphere start and end.

Usage

comp_ind_dj(d)

Arguments

Value

A vector of size sum(d) + 1.

Examples

# Example on (S^1)^3
d <- c(1, 1, 1)
comp_ind_dj(d = d)
comp_ind_dj(d = d) + 1

# Example on S^1 x S^2
d <- c(1, 2)
comp_ind_dj(d = d)
comp_ind_dj(d = d) + 1

Description

Computes the curvature matrix $\boldsymbol{R}(\boldsymbol{\kappa})$ of a product of von Mises–Fisher densities on the polysphere. This curvature is used in the rule-of-thumb selector bw_rot_polysph.

Usage

curv_vmf_polysph(kappa, d, log = FALSE)

Arguments

Value

A matrix of size c(length(r), length(r)).

Examples

# Curvature matrix
d <- 2:4
kappa <- 1:3
curv_vmf_polysph(kappa = kappa, d = d)
curv_vmf_polysph(kappa = kappa, d = d, log = TRUE)

# Equivalence on the sphere with DirStats::R_Psi_mixvmf
drop(curv_vmf_polysph(kappa = kappa[1], d = d[1]))
d[1]^2 * DirStats::R_Psi_mixvmf(q = d[1], mu = rbind(c(rep(0, d[1]), 1)),
                                kappa = kappa[1], p = 1)

Description

Computes the density of the uniform distribution on the polysphere.

Usage

d_unif_polysph(x, d, log = FALSE)

Arguments

Value

A vector of size nx with the evaluated density.

Examples

# Simple check of integration on S^1 x S^2
d <- c(1, 2)
x <- r_unif_polysph(n = 1e4, d = d)
mean(1 / d_unif_polysph(x = x, d = d)) / prod(rotasym::w_p(p = d + 1))

Description

Computes the density of the product of von Mises–Fisher densities on the polysphere.

Usage

d_vmf_polysph(x, d, mu, kappa, log = FALSE)

Arguments

Value

A vector of size nx with the evaluated density.

Examples

# Simple check of integration on S^1 x S^2
d <- c(1, 2)
mu <- c(0, 1, 0, 1, 0)
kappa <- c(1, 1)
x <- r_vmf_polysph(n = 1e4, d = d, mu = mu, kappa = kappa)
mean(1 / d_vmf_polysph(x = x, d = d, mu = mu, kappa = kappa)) /
  prod(rotasym::w_p(p = d + 1))

Description

Computation of the distance between points $\boldsymbol{x}$ and $\boldsymbol{y}$ on the polysphere $\mathcal{S}^{d_1} \times \cdots \times \mathcal{S}^{d_r}$:

$\sqrt{\sum_{j=1}^r d_{\mathcal{S}^{d_j}}(\boldsymbol{x}_j, \boldsymbol{y}_j)^2},$

where $d_{\mathcal{S}^{d_j}}(\boldsymbol{x}_j, \boldsymbol{y}_j)= \cos^{-1}(\boldsymbol{x}_j' \boldsymbol{y}_j)$.

Usage

dist_polysph(x, y, ind_dj, norm_x = FALSE, norm_y = FALSE, std = TRUE)

dist_polysph_cross(x, y, ind_dj, norm_x = FALSE, norm_y = FALSE,
  std = TRUE)

dist_polysph_matrix(x, ind_dj, norm_x = FALSE, norm_y = FALSE,
  std = TRUE)

Arguments

Value

• dist_polysph: a vector of size n with the distances between x and y.

• dist_polysph_matrix: a matrix of size c(n, n) with the pairwise distances of x.

• dist_polysph_cross: a matrix of distances of size c(n, m) with the cross distances between x and y.

Examples

# Example on S^2 x S^3 x S^1
d <- c(2, 3, 1)
ind_dj <- comp_ind_dj(d)
n <- 3
x <- r_unif_polysph(n = n, d = d)
y <- r_unif_polysph(n = n, d = d)

# Distances of x to y
dist_polysph(x = x, y = y, ind_dj = ind_dj, std = FALSE)
dist_polysph(x = x, y = y[1, , drop = FALSE], ind_dj = ind_dj, std = FALSE)

# Pairwise distance matrix of x
dist_polysph_matrix(x = x, ind_dj = ind_dj, std = FALSE)

# Cross distances between x and y
dist_polysph_cross(x = x, y = y, ind_dj = ind_dj, std = FALSE)

Description

Computes moments of kernels on $\mathcal{S}^{d_1} \times \cdots \times \mathcal{S}^{d_r}$ and efficiencies of kernels on $(\mathcal{S}^d)^r$.

Usage

eff_kern(d, r, k = 10, kernel, kernel_type = c("prod", "sph")[1],
  kernel_ref = "2", kernel_ref_type = c("prod", "sph")[2], ...)

b_d(kernel, d, k = 10, kernel_type = c("prod", "sph")[1], ...)

v_d(kernel, d, k = 10, kernel_type = c("prod", "sph")[1], ...)

Arguments

Value

• b_d: a vector with the first kernel moment on each hypersphere (common if kernel_type = "sph").

• v_d: a vector with the second kernel moment if kernel_type = "prod", or a scalar if kernel_type = "sph".

• eff_kern: a scalar with the kernel efficiency.

Examples

# Kernel moments
b_d(kernel = 2, d = c(2, 3), kernel_type = "prod")
v_d(kernel = 2, d = c(2, 3), kernel_type = "prod")
b_d(kernel = 2, d = c(2, 3), kernel_type = "sph")
v_d(kernel = 2, d = c(2, 3), kernel_type = "sph")

# Kernel efficiencies
eff_kern(d = 2, r = 1, kernel = "1")
eff_kern(d = 2, r = 1, kernel = "2")
eff_kern(d = 2, r = 1, k = 10, kernel = "3")

Description

Functions to perform density ridge estimation on the polysphere $\mathcal{S}^{d_1} \times \cdots \times \mathcal{S}^{d_r}$ through the Euler algorithm in standard, parallel, or block mode.

Usage

euler_ridge(x, X, d, h, h_euler = as.numeric(c()),
  weights = as.numeric(c()), wrt_unif = FALSE, normalized = TRUE,
  norm_x = FALSE, norm_X = FALSE, kernel = 1L, kernel_type = 1L,
  k = 10, N = 1000L, eps = 1e-05, keep_paths = FALSE,
  proj_alt = TRUE, fix_u1 = TRUE, sparse = FALSE, show_prog = TRUE,
  show_prog_j = FALSE)

parallel_euler_ridge(x, X, d, h, h_euler, N = 1000, eps = 1e-05,
  keep_paths = FALSE, cores = 1, ...)

block_euler_ridge(x, X, d, h, h_euler, ind_blocks, N = 1000, eps = 1e-05,
  keep_paths = FALSE, cores = 1, ...)

Arguments

Details

euler_ridge is the main function to perform density ridge estimation through the Euler algorithm from the starting values x to initiate the ridge path. The function euler_ridge_parallel parallelizes on the starting values x. The function euler_ridge_block runs the Euler algorithm marginally in blocks of hyperspheres, instead of jointly in the whole polysphere. This function requires that all the dimensions are the same.

Value

The three functions return a list with the following fields:

Examples

## Test on S^2 with a small circle trend

# Sample
r <- 1
d <- 2
n <- 50
ind_dj <- comp_ind_dj(d = d)
set.seed(987204452)
X <- r_path_s2r(n = n, r = r, spiral = FALSE, Theta = cbind(c(1, 0, 0)),
                sigma = 0.35)[, , 1]
col_X_alp <- viridis::viridis(n, alpha = 0.25)
col_X <- viridis::viridis(n)

# Euler
h_rid <- 0.5
h_eu <- h_rid^2
N <- 30
eps <- 1e-6
Y <- euler_ridge(x = X, X = X, d = d, h = h_rid, h_euler = h_eu,
                 N = N, eps = eps, keep_paths = TRUE)
Y

# Visualization
i <- N # Between 1 and N
sc3 <- scatterplot3d::scatterplot3d(Y$paths[, , 1], color = col_X_alp,
                                    pch = 19, xlim = c(-1, 1),
                                    ylim = c(-1, 1), zlim = c(-1, 1),
                                    xlab = "x", ylab = "y", zlab = "z")
sc3$points3d(rbind(Y$paths[, , i]), col = col_X, pch = 16, cex = 0.75)
for (k in seq_len(nrow(Y$paths))) {

  sc3$points3d(t(Y$paths[k, , ]), col = col_X_alp[k], type = "l")

}

Description

Computes the gradient $\mathsf{D}\hat{f}(\boldsymbol{x};\boldsymbol{h})$ and Hessian matrix $\mathsf{H}\hat{f}(\boldsymbol{x};\boldsymbol{h})$ of the kernel density estimator $\hat{f}(\boldsymbol{x};\boldsymbol{h})$ on the polysphere $\mathcal{S}^{d_1} \times \cdots \times \mathcal{S}^{d_r}$.

Usage

grad_hess_kde_polysph(x, X, d, h, weights = as.numeric(c()),
  projected = TRUE, proj_alt = TRUE, norm_grad_hess = FALSE,
  log = FALSE, wrt_unif = FALSE, normalized = TRUE, norm_x = FALSE,
  norm_X = FALSE, kernel = 1L, kernel_type = 1L, k = 10)

Arguments

Value

A list with the following components:

Examples

# Simple check on (S^1)^2
n <- 3
d <- c(1, 1)
mu <- c(0, 1, 0, 1)
kappa <- c(5, 5)
h <- c(0.2, 0.2)
X <- r_vmf_polysph(n = n, d = d, mu = mu, kappa = kappa)
grh <- grad_hess_kde_polysph(x = X, X = X, d = d, h = h)
str(grh)
grh

Description

Permutation tests for the equality of distributions of two or $k$ samples of data on $\mathcal{S}^{d_1} \times \cdots \times \mathcal{S}^{d_r}$. The Jensen–Shannon distance is used to construct a test statistic measuring the discrepancy between the $k$ kernel density estimators. Tests based on the mean and scatter matrices are also available, but for only two samples ($k=2$).

Usage

hom_test_polysph(X, d, labels, type = c("jsd", "mean", "scatter", "hd")[1],
  h = NULL, kernel = 1, kernel_type = 1, k = 10, B = 1000,
  M = 10000, plot_boot = FALSE, seed_jsd = NULL, cv_jsd = TRUE)

Arguments

Details

Only type = "jsd" is able to deal with $k > 2$.

The "jsd" statistic is the Jensen–Shannon divergence. This statistic is bounded in $[0, 1]$. The "mean" statistic measures the maximum (chordal) distance between the estimated group means. This statistic is bounded in $[0, 1]$. The "scatter" statistic measures the maximum affine invariant Riemannian metric between the estimated scatter matrices. The "hd" statistic computes a monotonic transformation of the Hellinger distance, which is the Bhattacharyya divergence (or coefficient).

Value

An object of class "htest" with the following fields:

Examples

## Two-sample case

# H0 holds
n <- c(50, 100)
X1 <- rotasym::r_vMF(n = n[1], mu = c(0, 0, 1), kappa = 1)
X2 <- rotasym::r_vMF(n = n[2], mu = c(0, 0, 1), kappa = 1)
hom_test_polysph(X = rbind(X1, X2), labels = rep(1:2, times = n),
                 d = 2, type = "jsd", h = 0.5)

# H0 does not hold
X2 <- rotasym::r_vMF(n = n[2], mu = c(0, 1, 0), kappa = 2)
hom_test_polysph(X = rbind(X1, X2), labels = rep(1:2, times = n),
                 d = 2, type = "jsd", h = 0.5)

## k-sample case

# H0 holds
n <- c(50, 100, 50)
X1 <- rotasym::r_vMF(n = n[1], mu = c(0, 0, 1), kappa = 1)
X2 <- rotasym::r_vMF(n = n[2], mu = c(0, 0, 1), kappa = 1)
X3 <- rotasym::r_vMF(n = n[3], mu = c(0, 0, 1), kappa = 1)
hom_test_polysph(X = rbind(X1, X2, X3), labels = rep(1:3, times = n),
                 d = 2, type = "jsd", h = 0.5)

# H0 does not hold
X3 <- rotasym::r_vMF(n = n[3], mu = c(0, 1, 0), kappa = 2)
hom_test_polysph(X = rbind(X1, X2, X3), labels = rep(1:3, times = n),
                 d = 2, type = "jsd", h = 0.5)

Description

Indexing of an unordered collection of points defining the estimated density ridge curve. The indexing is done by a multidimensional scaling map to the real line, while the smoothing is done by local polynomial regression for polyspherical-on-scalar regression.

Usage

index_ridge(endpoints, X, d, l_index = 1000, f_index = 2,
  probs_scores = seq(0, 1, l = 101), verbose = FALSE, type_bwd = c("1se",
  "min")[1], p = 0, ...)

Arguments

Details

Indexing is designed to work with collection of ridge points that admit a linear ordering, so that mapping to the real line by multidimensional scaling is adequate. The indexing will not work properly if the ridge points define a closed curve.

Value

A list with the following fields:

Examples

## Test on (S^1)^2

# Sample
set.seed(132121)
r <- 2
d <- rep(1, r)
n <- 200
ind_dj <- comp_ind_dj(d = d)
Th <- matrix(runif(n = n * (r - 1), min = -pi / 2, max = pi / 2),
             nrow = n, ncol = r - 1)
Th[, r - 1] <- sort(Th[, r - 1])
Th <- cbind(Th, sdetorus::toPiInt(
  pi + Th[, r - 1] + runif(n = n, min = -pi / 4, max = pi / 4)))
X <- angles_to_torus(Th)
col_X_alp <- viridis::viridis(n, alpha = 0.25)
col_X <- viridis::viridis(n)

# Euler
h_rid <- rep(0.75, r)
h_eu <- h_rid^2
N <- 200
eps <- 1e-6
Y <- euler_ridge(x = X, X = X, d = d, h = h_rid, h_euler = h_eu,
                 N = N, eps = eps, keep_paths = TRUE)

# Visualization
i <- N # Between 1 and N
plot(rbind(torus_to_angles(Y$paths[, , 1])), col = col_X_alp, pch = 19,
     axes = FALSE, xlim = c(-pi, pi), ylim = c(-pi, pi),
     xlab = "", ylab = "")
points(rbind(torus_to_angles(Y$paths[, , i])), col = col_X, pch = 16,
       cex = 0.75)
sdetorus::torusAxis(1:2)
for (k in seq_len(nrow(Y$paths))) {

  xy <- torus_to_angles(t(Y$paths[k, , ]))
  sdetorus::linesTorus(x = xy[, 1], y = xy[, 2], col = col_X_alp[k])

}

# SIER
ind_rid <- index_ridge(endpoints = Y$ridge_y, X = X, d = d,
                       probs_scores = seq(0, 1, l = 50))
xy <- torus_to_angles(ind_rid$ridge_grid)
sdetorus::linesTorus(x = xy[, 1], y = xy[, 2], col = 2, lwd = 2)
points(torus_to_angles(ind_rid$ridge_fun(quantile(ind_rid$scores_grid))),
       col = 4, pch = 19)

# Scores
plot(density(ind_rid$scores_X), type = "l", xlab = "Scores",
     ylab = "Density", main = "Scores density")
for (i in 1:n) rug(ind_rid$scores_X[i], col = col_X[i])

Description

Creates a sequence of points on the polysphere linearly interpolating between two points extrinsically.

Usage

interp_polysph(x, y, ind_dj, N = 10)

Arguments

Value

A matrix of size c(N, sum(d) + r) with the interpolation.

Examples

interp_polysph(x = c(1, 0), y = c(0, 1), ind_dj = comp_ind_dj(d = 1))

Description

Computes the kernel density estimator for data on the polysphere $\mathcal{S}^{d_1} \times \cdots \times \mathcal{S}^{d_r}$. Given a sample $\boldsymbol{X}_1,\ldots,\boldsymbol{X}_n$, this estimator is

$\hat{f}(\boldsymbol{x};\boldsymbol{h})=\sum_{i=1}^n L_{\boldsymbol{h}}(\boldsymbol{x},\boldsymbol{X}_i)$

for a kernel $L$ and a vector of bandwidths $\boldsymbol{h}$.

Usage

kde_polysph(x, X, d, h, weights = as.numeric(c()), log = FALSE,
  wrt_unif = FALSE, normalized = TRUE, intrinsic = FALSE,
  norm_x = FALSE, norm_X = FALSE, kernel = 1L, kernel_type = 1L,
  k = 10)

Arguments

Value

A column vector of size c(nx, 1) with the evaluation of kernel density estimator.

Examples

# Simple check on S^1 x S^2
n <- 1e3
d <- c(1, 2)
mu <- c(0, 1, 0, 0, 1)
kappa <- c(5, 5)
h <- c(0.2, 0.2)
X <- r_vmf_polysph(n = n, d = d, mu = mu, kappa = kappa)
kde_polysph(x = rbind(mu), X = X, d = d, h = h)
d_vmf_polysph(x = rbind(mu), d = d, mu = mu, kappa = kappa)

Description

An isotropic kernel $L$ on $\mathcal{S}^d$ and its normalizing constant are such that $\int_{\mathcal{S}^d} c(h, d, L) L\left(\frac{1 - \boldsymbol{x}'\boldsymbol{y}}{h^2}\right) \,\mathrm{d}\boldsymbol{x} = 1$ (extrinsic-chordal distance) or $\int_{\mathcal{S}^d} c(h, d, L) L\left(\frac{\cos^{-1}(\boldsymbol{x}'\boldsymbol{y})^2}{2h^2}\right) \,\mathrm{d}\boldsymbol{x} = 1$ (intrinsic distance).

Usage

L(t, kernel = "1", squared = FALSE, deriv = 0, k = 10,
  inc_sfp = TRUE)

c_kern(h, d, kernel = "1", kernel_type = "1", k = 10, log = FALSE,
  inc_sfp = TRUE, intrinsic = FALSE)

grad_L(x, y, h, kernel = 1, k = 10)

hess_L(x, y, h, kernel = 1, k = 10)

Arguments

Details

The gradient and Hessian are computed for the functions $\boldsymbol{x} \mapsto L\left(\frac{1 - \boldsymbol{x}'\boldsymbol{y}}{h^2}\right)$.

Value

• L: a vector with the kernel evaluated at t.

• grad_L: a vector with the gradient evaluated at x.

• hess_L: a matrix with the Hessian evaluated at x.

Examples

# Constants in terms of h
h_grid <- seq(0.01, 4, l = 100)
r <- 2
d <- 2
dr <- rep(d, r)
c_vmf <- sapply(h_grid, function(hi)
  log(c_kern(h = rep(hi, r), d = dr, kernel = 1, kernel_type = 2)))
c_epa <- sapply(h_grid, function(hi)
  log(c_kern(h = rep(hi, r), d = dr, kernel = 2, kernel_type = 2)))
c_sfp <- sapply(h_grid, function(hi)
  log(c_kern(h = rep(hi, r), d = dr, kernel = 3, k = 1, kernel_type = 2)))
plot(h_grid, c_epa, type = "l", ylab = "Constant", xlab = "h", col = 2)
lines(h_grid, c_sfp, col = 3)
lines(h_grid, c_vmf, col = 1)
abline(v = sqrt(2), lty = 2, col = 2)

# Kernel and its derivatives
h <- 0.5
x <- c(sqrt(2), -sqrt(2), 0) / 2
y <- c(-sqrt(2), sqrt(3), sqrt(3)) / 3
L(t = (1 - sum(x * y)) / h^2)
grad_L(x = x, y = y, h = h)
hess_L(x = x, y = y, h = h)

Description

Computes a local constant (Nadaraya–Watson) or local linear estimator with polyspherical response and scalar predictor.

Usage

kre_polysph(x, X, Y, d, h, p = 0)

Arguments

Value

A vector of size nx with the estimated regression curve evaluated at x.

Examples

x_grid <- seq(-0.25, 1.25, l = 200)
n <- 50
X <- seq(0, 1, l = n)
Y <- r_path_s2r(n = n, r = 1, sigma = 0.1, spiral = TRUE)[, , 1]
h0 <- bw_cv_kre_polysph(X = X, Y = Y, d = 2, p = 0, plot_cv = FALSE)$h_1se
sc3 <- scatterplot3d::scatterplot3d(Y, pch = 16, xlim = c(-1, 1),
                                    ylim = c(-1, 1), zlim = c(-1, 1),
                                    xlab = "", ylab = "", zlab = "")
sc3$points3d(kre_polysph(x = x_grid, X = X, Y = Y, d = 2, h = h0, p = 0),
             pch = 16, type = "l", col = 2, lwd = 2)
sc3$points3d(kre_polysph(x = x_grid, X = X, Y = Y, d = 2, h = h0, p = 1),
             pch = 16, type = "l", col = 3, lwd = 2)

Description

Computes the logarithm of the cross-validated kernel density estimator: $\log \hat{f}_{-i}(\boldsymbol{X}_i;\boldsymbol{h})$, $i = 1, \ldots, n.$

Usage

log_cv_kde_polysph(X, d, h, weights = as.numeric(c()), wrt_unif = FALSE,
  normalized = TRUE, intrinsic = FALSE, norm_X = FALSE, kernel = 1L,
  kernel_type = 1L, k = 10)

Arguments

Value

A column vector of size c(n, 1) with the evaluation of the logarithm of the cross-validated kernel density estimator.

Examples

# Simple check on S^1 x S^2
n <- 5
d <- c(1, 2)
h <- c(0.2, 0.2)
X <- r_unif_polysph(n = n, d = d)
log_cv_kde_polysph(X = X, d = d, h = h)
kde_polysph(x = X[1, , drop = FALSE], X = X[-1, ], d = d, h = h, log = TRUE)

Description

Computation of the polylogarithm $\mathrm{Li}_s(-e^\mu)$.

Usage

polylog_minus_exp_mu(mu, s, upper = Inf, ...)

Arguments

Details

If s is an integer, 1/2, 3/2, or 5/2, then routines from the GSL library to compute Fermi–Dirac integrals are called. Otherwise, numerical integration is used.

Value

A vector of size length(mu) or length(s) with the values of the polylogarithm.

Examples

polylog_minus_exp_mu(mu = 1:5, s = 1)
polylog_minus_exp_mu(mu = 1, s = 1:5)
polylog_minus_exp_mu(mu = 1:5, s = 1:5)

Description

Computes the projected gradient $\mathsf{D}_{(p-1)}\hat{f}(\boldsymbol{x};\boldsymbol{h})$ of the kernel density estimator $\hat{f}(\boldsymbol{x};\boldsymbol{h})$ on the polysphere $\mathcal{S}^{d_1} \times \cdots \times \mathcal{S}^{d_r}$, where $p=\sum_{j=1}^r d_j+r$ is the dimension of the ambient space.

Usage

proj_grad_kde_polysph(x, X, d, h, weights = as.numeric(c()),
  wrt_unif = FALSE, normalized = TRUE, norm_x = FALSE, norm_X = FALSE,
  kernel = 1L, kernel_type = 1L, k = 10, proj_alt = TRUE,
  fix_u1 = TRUE, sparse = FALSE)

Arguments

Value

A list with the following components:

Examples

# Simple check on (S^1)^2
n <- 3
d <- c(1, 1)
mu <- c(0, 1, 0, 1)
kappa <- c(5, 5)
h <- c(0.2, 0.2)
X <- r_vmf_polysph(n = n, d = d, mu = mu, kappa = kappa)
proj_grad_kde_polysph(x = X, X = X, d = d, h = h)

Description

Projects points on $\mathbb{R}^{d_1 + \cdots + d_r + r}$ onto the polysphere $\mathcal{S}^{d_1} \times \cdots \times \mathcal{S}^{d_r}$ by normalizing each block of $d_j$ coordinates.

Usage

proj_polysph(x, ind_dj)

Arguments

Value

A matrix of size c(n, sum(d) + r) with the projected points.

Examples

# Example on (S^1)^2
d <- c(1, 1)
x <- rbind(c(2, 0, 1, 1))
proj_polysph(x, ind_dj = comp_ind_dj(d))

Description

Simulation from the angular density function of an isotropic kernel on the hypersphere $\mathcal{S}^d$.

Usage

r_g_kern(n, d, h, kernel = "1", k = 10)

Arguments

Value

A vector of size n with the sample.

Examples

hist(r_g_kern(n = 1e3, d = 2, h = 1, kernel = "1"), breaks = 30,
     probability = TRUE, main = "", xlim = c(-1, 1))
hist(r_g_kern(n = 1e3, d = 2, h = 1, kernel = "2"), breaks = 30,
     probability = TRUE, main = "", xlim = c(-1, 1))
hist(r_g_kern(n = 1e3, d = 2, h = 1, kernel = "3"), breaks = 30,
     probability = TRUE, main = "", xlim = c(-1, 1))

Description

Simulates from the distribution defined by a polyspherical kernel density estimator on $\mathcal{S}^{d_1} \times \ldots \times \mathcal{S}^{d_r}$.

Usage

r_kde_polysph(n, X, d, h, kernel = 1, kernel_type = 1, k = 10,
  norm_X = FALSE)

Arguments

Details

The function uses r_kern_polysph to sample from the considered kernel.

Value

A matrix of size c(n, sum(d) + r) with the sample.

Examples

# Simulated data on (S^1)^2
n <- 50
samp <- r_path_s1r(n = n, r = 2, k = c(1, 2), angles = TRUE)
plot(samp, xlim = c(-pi, pi), ylim = c(-pi, pi), col = rainbow(n),
     axes = FALSE, xlab = "", ylab = "", pch = 16, cex = 0.75)
points(torus_to_angles(r_kde_polysph(n = 10 * n, X = angles_to_torus(samp),
                                     d = c(1, 1), h = c(0.1, 0.1))),
       col = "black", pch = 16, cex = 0.2)
sdetorus::torusAxis()

# Simulated data on S^2
n <- 50
samp <- r_path_s2r(n = n, r = 1, sigma = 0.1, kappa = 5,
                   spiral = TRUE)[, , 1]
sc3d <- scatterplot3d::scatterplot3d(
  samp, xlim = c(-1, 1), ylim = c(-1, 1), zlim = c(-1, 1),
  xlab = "", ylab = "", zlab = "", color = rainbow(n), pch = 16
)
xyz <- r_kde_polysph(n = 10 * n, X = samp, d = 2, h = 0.1)
sc3d$points3d(xyz[, 1], xyz[, 2], xyz[, 3], col = "black", pch = 16,
              cex = 0.2)

Description

Simulates from the distribution defined by a kernel on the polysphere $\mathcal{S}^{d_1} \times \cdots \times \mathcal{S}^{d_r}$.

Usage

r_kern_polysph(n, d, mu, h, kernel = 1, kernel_type = 1, k = 10,
  norm_mu = FALSE)

Arguments

Details

Simulation for non-von Mises–Fisher spherically symmetric kernels is done by acceptance-rejection from a von Mises–Fisher proposal distribution.

Value

A matrix of size c(n, sum(d) + r) with the sample.

Examples

# Simulate kernels in (S^1)^2
n <- 1e3
h <- c(1, 1)
d <- c(1, 1)
mu <- rep(DirStats::to_cir(pi), 2)
samp_ker <- function(kernel, kernel_type, col, main) {
  data <- r_kern_polysph(n = n, d = d, mu = mu, h = h, kernel = kernel,
                         kernel_type = kernel_type)
  ang <- cbind(DirStats::to_rad(data[, 1:2]),
               DirStats::to_rad(data[, 3:4]))
  plot(ang, xlim = c(0, 2 * pi), ylim = c(0, 2 * pi), pch = 16, cex = 0.25,
       col = col, xlab = expression(Theta[1]), ylab = expression(Theta[2]),
       main = main)
}
old_par <- par(mfcol = c(2, 3))
samp_ker(kernel = 2, kernel_type = 2, col = 1, main = "Epa sph. symmetric")
samp_ker(kernel = 2, kernel_type = 1, col = 2, main = "Epa product")
samp_ker(kernel = 3, kernel_type = 2, col = 1, main = "Sfp sph. symmetric")
samp_ker(kernel = 3, kernel_type = 1, col = 2, main = "Sfp product")
samp_ker(kernel = 1, kernel_type = 2, col = 1, main = "vMF sph. symmetric")
samp_ker(kernel = 1, kernel_type = 1, col = 2, main = "vMF product")
par(old_par)

# Simulate kernels in (S^2)^2
n <- 1e3
h <- c(0.2, 0.6)
d <- c(2, 2)
mu <- c(c(0, 0, 1), c(0, -1, 0))
samp_ker <- function(kernel, kernel_type, main) {
  data <- r_kern_polysph(n = n, d = d, mu = mu, h = h, kernel = kernel,
                         kernel_type = kernel_type)
  scatterplot3d::scatterplot3d(rbind(data[, 1:3], data[, 4:6]),
                               xlim = c(-1, 1), ylim = c(-1, 1),
                               zlim = c(-1, 1), pch = 16, xlab = "",
                               ylab = "", zlab = "", cex.symbols = 0.5,
       color = rep(viridis::viridis(n)[rank(data[, 3])], 2), main = main)
}
old_par <- par(mfcol = c(2, 3))
samp_ker(kernel = 2, kernel_type = 2, main = "Epa sph. symmetric")
samp_ker(kernel = 2, kernel_type = 1, main = "Epa product")
samp_ker(kernel = 3, kernel_type = 2, main = "Sfp sph. symmetric")
samp_ker(kernel = 3, kernel_type = 1, main = "Sfp product")
samp_ker(kernel = 1, kernel_type = 2, main = "vMF sph. symmetric")
samp_ker(kernel = 1, kernel_type = 1, main = "vMF product")
par(old_par)

# Plot simulated data
n <- 1e3
h <- c(1, 1)
d <- c(2, 2)
samp_ker <- function(kernel, kernel_type, col, main) {
  X <- r_kern_polysph(n = n, d = d, mu = mu, h = h, kernel = kernel,
                      kernel_type = kernel_type)
  S <- cbind((1 - X[, 1:3] %*% mu[1:3]) / h[1]^2,
             (1 - X[, 4:6] %*% mu[4:6]) / h[2]^2)
  plot(S, xlim = c(0, 2 / h[1]^2), ylim = c(0, 2 / h[2]^2), pch = 16,
       cex = 0.25, col = col, xlab = expression(t[1]),
       ylab = expression(t[2]), main = main)
  t_grid <- seq(0, 2 / min(h)^2, l = 100)
  gr <- as.matrix(expand.grid(t_grid, t_grid))
  if (kernel_type == "1") {

    dens <- prod(c_kern(h = h, d = d, kernel = kernel, kernel_type = 1)) *
      L(gr[, 1], kernel = kernel) * L(gr[, 2], kernel = kernel)

  } else if (kernel_type == "2") {

    dens <- c_kern(h = h, d = d, kernel = kernel, kernel_type = 2) *
      L(gr[, 1] + gr[, 2], kernel = kernel)

  }
  dens <- matrix(dens, nrow = length(t_grid), ncol = length(t_grid))
  contour(t_grid, t_grid, dens, add = TRUE, col = col,
          levels = seq(0, 0.2, l = 41))

}
old_par <- par(mfcol = c(2, 3))
samp_ker(kernel = 2, kernel_type = 2, col = 1, main = "Epa sph. symmetric")
samp_ker(kernel = 2, kernel_type = 1, col = 2, main = "Epa product")
samp_ker(kernel = 3, kernel_type = 2, col = 1, main = "Sfp sph. symmetric")
samp_ker(kernel = 3, kernel_type = 1, col = 2, main = "Sfp product")
samp_ker(kernel = 1, kernel_type = 2, col = 1, main = "vMF sph. symmetric")
samp_ker(kernel = 1, kernel_type = 1, col = 2, main = "vMF product")
par(old_par)