euler_ridge: Euler algorithms for polyspherical density ridge estimation

In polykde: Polyspherical Kernel Density Estimation

Euler algorithms for polyspherical density ridge estimation

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

x	a matrix of size c(nx, sum(d) + r) with the starting points for the Euler algorithm.
X	a matrix of size c(n, sum(d) + r) with the sample.
d	vector of size r with dimensions.
h	vector of size r with bandwidths.
h_euler	vector of size r with the advance steps in the Euler method. Set internally as h if not provided.
weights	weights for each observation. If provided, a vector of size n with the weights for multiplying each kernel. If not provided, set internally to rep(1 / n, n), which gives the standard estimator.
wrt_unif	flag to return a density with respect to the uniform measure. If FALSE (default), the density is with respect to the Lebesgue measure.
normalized	flag to compute the normalizing constant of the kernel and include it in the kernel density estimator. Defaults to TRUE.
norm_x, norm_X	ensure a normalization of the data? Defaults to FALSE.
kernel	kernel employed: 1 for von Mises–Fisher (default); 2 for Epanechnikov; 3 for softplus.
kernel_type	type of kernel employed: 1 for product kernel (default); 2 for spherically symmetric kernel.
k	softplus kernel parameter. Defaults to 10.0.
N	maximum number of Euler iterations. Defaults to 1e3.
eps	convergence tolerance. Defaults to 1e-5.
keep_paths	keep the Euler paths to the ridge? Defaults to FALSE.
proj_alt	alternative projection. Defaults to TRUE.
fix_u1	ensure the u_1 vector is different from x? Prevents the Euler algorithm to "surf the ridge". Defaults to TRUE.
sparse	use a sparse eigendecomposition of the Hessian? Defaults to FALSE.
show_prog	display a progress bar for x? Defaults to TRUE.
show_prog_j	display a progress bar for N? Defaults to FALSE.
cores	cores to use. Defaults to 1.
...	further arguments passed to euler_ridge.
ind_blocks	indexes of the blocks, a vector or length r.

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 spheres, 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:

ridge_y	a matrix of size c(nx, sum(d) + r) with the end points of Euler algorithm defining the estimated ridge.
lamb_norm_y	a matrix of size c(nx, sum(d) + r) with the Hessian eigenvalues (largest to smallest) evaluated at end points.
log_dens_y	a column vector of size c(nx, 1) with the logarithm of the density at end points.
paths	an array of size c(nx, sum(d) + r, N + 1) containing the Euler paths.
start_x	a matrix of size c(nx, sum(d) + r) with the starting points for the Euler algorithm.
iter	a column vector of size c(nx, 1) counting the iterations required for each point.
conv	a column vector of size c(nx, 1) with convergence flags.
d	vector d.
h	bandwidth used for the kernel density estimator.
error	a column vector of size c(nx, 1) indicating if errors were found for each path.

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")

}

polykde documentation built on April 16, 2025, 1:11 a.m.