tang-norm-decomp: Distributions based on the tangent-normal decomposition

Description

Density and simulation of a distribution on $S^{p-1}:=\{\mathbf{x}\in R^p:||\mathbf{x}||=1\}$, $p\ge 2$, obtained by the tangent-normal decomposition. The tangent-normal decomposition of the random vector $\mathbf{X}\in S^{p-1}$ is $$V\boldsymbol{\theta} + \sqrt{1 - V^2}\boldsymbol{\Gamma}_{\boldsymbol{\theta}}\mathbf{U}$$ where $V := \mathbf{X}'\boldsymbol{\theta}$ is a random variable in $[-1, 1]$ (the cosines of $\mathbf{X}$) and $\mathbf{U} := \boldsymbol{\Gamma}_{\boldsymbol{\theta}}\mathbf{X}/ ||\boldsymbol{\Gamma}_{\boldsymbol{\theta}}\mathbf{X}||$ is a random vector in $S^{p-2}$ (the multivariate signs of $\mathbf{X}$) and $\boldsymbol{\Gamma}_{\boldsymbol{\theta}}$ is the $p\times(p-1)$ matrix computed by Gamma_theta.

The tangent-normal decomposition can be employed for constructing distributions for $\mathbf{X}$ that arise for certain choices of $V$ and $\mathbf{U}$. If $V$ and $\mathbf{U}$ are independent, then simulation from $\mathbf{X}$ is straightforward using the tangent-normal decomposition. Also, the density of $\mathbf{X}$ at $\mathbf{x}\in S^{p-1}$, $f_\mathbf{X}(\mathbf{x})$, is readily computed as $$f_\mathbf{X}(\mathbf{x})= \omega_{p-1}c_g g(t)(1-t^2)^{(p-3)/2}f_\mathbf{U}(\mathbf{u})$$ where $t:=\mathbf{x}'\boldsymbol{\theta}$, $\mathbf{u}:=\boldsymbol{\Gamma}_{\boldsymbol{\theta}}\mathbf{x}/ ||\boldsymbol{\Gamma}_{\boldsymbol{\theta}}\mathbf{x}||$, $f_\mathbf{U}$ is the density of $\mathbf{U}$, and $f_V(v) := \omega_{p-1} c_g g(v) (1 - v^2)^{(p-3)/2}$ is the density of $V$ for an angular function $g$ with normalizing constant $c_g$. $\omega_{p-1}$ is the surface area of $S^{p-2}$.

Usage

d_tang_norm(x, theta, g_scaled, d_V, d_U, log = FALSE)
r_tang_norm(n, theta, r_U, r_V)

Value

Depending on the function:

d_tang_norm: a vector of length nx or 1 with the evaluated density at x.
r_tang_norm: a matrix of size c(n, p) with the random sample.

Arguments

x: locations in $S^{p-1}$ to evaluate the density. Either a matrix of size c(nx, p) or a vector of length p. Normalized internally if required (with a warning message).
theta: a unit norm vector of size p giving the axis of rotational symmetry.
g_scaled: the scaled angular density $c_g g$. In the form
g_scaled <- function(t, log = TRUE) {...}. See examples.
d_V: the density $f_V$. In the form d_V <- function(v, log = TRUE) {...}. See examples.
d_U: the density $f_\mathbf{U}$. In the form d_U <- function(u, log = TRUE) {...}. See examples.
log: flag to indicate if the logarithm of the density (or the normalizing constant) is to be computed.
n: sample size, a positive integer.
r_U: a function for simulating $\mathbf{U}$. Its first argument must be the sample size. See examples.
r_V: a function for simulating $V$. Its first argument must be the sample size. See examples.

Author

Eduardo García-Portugués, Davy Paindaveine, and Thomas Verdebout.

Details

Either g_scaled or d_V can be supplied to d_tang_norm (the rest of the arguments are compulsory). One possible choice for g_scaled is g_vMF with scaled = TRUE. Another possible choice is the angular function $g(t) = 1 - t^2$, normalized by its normalizing constant $c_g = (\Gamma(p/2) p) / (2\pi^{p/2} (p - 1))$ (see examples). This angular function makes $V^2$ to be distributed as a $\mathrm{Beta}(1/2,(p+1)/2)$.

The normalizing constants and densities are computed through log-scales for numerical accuracy.

References

García-Portugués, E., Paindaveine, D., Verdebout, T. (2020) On optimal tests for rotational symmetry against new classes of hyperspherical distributions. Journal of the American Statistical Association, 115(532):1873--1887. tools:::Rd_expr_doi("10.1080/01621459.2019.1665527")

Examples

Run this code

## Simulation and density evaluation for p = 2

# Parameters
n <- 1e3
p <- 2
theta <- c(rep(0, p - 1), 1)
mu <- c(rep(0, p - 2), 1)
kappa_V <- 2
kappa_U <- 0.1

# The vMF scaled angular function
g_scaled <- function(t, log) {
  g_vMF(t, p = p - 1, kappa = kappa_V, scaled = TRUE, log = log)
}

# Cosine density for the vMF distribution
d_V <- function(v, log) {
 log_dens <- g_scaled(v, log = log) + (p - 3)/2 * log(1 - v^2)
 switch(log + 1, exp(log_dens), log_dens)
}

# Multivariate signs density based on a vMF
d_U <- function(x, log) d_vMF(x = x, mu = mu, kappa = kappa_U, log = log)

# Simulation functions
r_V <- function(n) r_g_vMF(n = n, p = p, kappa = kappa_V)
r_U <- function(n) r_vMF(n = n, mu = mu, kappa = kappa_U)

# Sample and color according to density
x <- r_tang_norm(n = n, theta = theta, r_V = r_V, r_U = r_U)
r <- runif(n, 0.95, 1.05) # Radius perturbation to improve visualization
col <- viridisLite::viridis(n)
dens <- d_tang_norm(x = x, theta = theta, g_scaled = g_scaled, d_U = d_U)
# dens <- d_tang_norm(x = x, theta = theta, d_V = d_V, d_U = d_U) # The same
plot(r * x, pch = 16, col = col[rank(dens)])

## Simulation and density evaluation for p = 3

# Parameters
p <- 3
n <- 5e3
theta <- c(rep(0, p - 1), 1)
mu <- c(rep(0, p - 2), 1)
kappa_V <- 2
kappa_U <- 2

# Sample and color according to density
x <- r_tang_norm(n = n, theta = theta, r_V = r_V, r_U = r_U)
col <- viridisLite::viridis(n)
dens <- d_tang_norm(x = x, theta = theta, g_scaled = g_scaled, d_U = d_U)
if (requireNamespace("rgl")) {
  rgl::plot3d(x, col = col[rank(dens)], size = 5)
}

## A non-vMF angular function: g(t) = 1 - t^2. It is sssociated to the
## Beta(1/2, (p + 1)/2) distribution.

# Scaled angular function
g_scaled <- function(t, log) {
  log_c_g <- lgamma(0.5 * p) + log(0.5 * p / (p - 1)) - 0.5 * p * log(pi)
  log_g <- log_c_g + log(1 - t^2)
  switch(log + 1, exp(log_g), log_g)
}

# Cosine density
d_V <- function(v, log) {
  log_dens <- w_p(p = p - 1, log = TRUE) + g_scaled(t = v, log = TRUE) +
    (0.5 * (p - 3)) * log(1 - v^2)
  switch(log + 1, exp(log_dens), log_dens)
}

# Simulation
r_V <- function(n) {
  sample(x = c(-1, 1), size = n, replace = TRUE) *
    sqrt(rbeta(n = n, shape1 = 0.5, shape2 = 0.5 * (p + 1)))
}

# Sample and color according to density
r_U <- function(n) r_unif_sphere(n = n, p = p - 1)
x <- r_tang_norm(n = n, theta = theta, r_V = r_V, r_U = r_U)
col <- viridisLite::viridis(n)
dens <- d_tang_norm(x = x, theta = theta, d_V = d_V, d_U = d_unif_sphere)
# dens <- d_tang_norm(x = x, theta = theta, g_scaled = g_scaled,
#                     d_U = d_unif_sphere) # The same
if (requireNamespace("rgl")) {
  rgl::plot3d(x, col = col[rank(dens)], size = 5)
}

Run the code above in your browser using DataLab