polykde (version 1.1.7)

euler_ridge: 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, ...)

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.

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.

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

}