# R/ACG.R In rotasym: Tests for Rotational Symmetry on the Hypersphere

#### Documented in c_ACGd_ACGr_ACG

#' @title Angular central Gaussian distribution
#'
#' @description Density and simulation of the Angular Central Gaussian (ACG)
#' distribution on
#' \eqn{S^{p-1}:=\{\mathbf{x}\in R^p:||\mathbf{x}||=1\}}{
#' S^{p-1} := \{x \in R^p : ||x|| = 1\}}, \eqn{p\ge 1}. The density at
#' \eqn{\mathbf{x} \in S^{p-1}}{x \in S^{p-1}}, \eqn{p\ge 2}, is given by
#' \deqn{c^{\mathrm{ACG}}_{p,\boldsymbol{\Lambda}}
#' (\mathbf{x}' \boldsymbol{\Lambda}^{-1} \mathbf{x})^{-p/2}
#' 1 / (\omega_p |\boldsymbol{\Lambda}|^{1/2})}{
#' c^{ACG}_{p,\Lambda} (x' \Lambda^{-1} x)^{-p/2} with
#' c^{ACG}_{p,\Lambda} := 1 / (\omega_p |\Lambda|^{1/2}).}
#' where \eqn{\boldsymbol{\Lambda}}{\Lambda} is the shape matrix, a
#' \eqn{p\times p}{p x p} symmetric and positive definite matrix, and
#' \eqn{\omega_p} is the surface area of \eqn{S^{p-1}}.
#'
#' @inheritParams unif
#' @param Lambda the shape matrix \eqn{\boldsymbol{\Lambda}}{\Lambda} of the
#' ACG. A symmetric and positive definite matrix of size \code{c(p, p)}.
#' @return
#' Depending on the function:
#' \itemize{
#' \item \code{d_ACG}: a vector of length \code{nx} or \code{1} with the
#' evaluated density at \code{x}.
#' \item \code{r_ACG}: a matrix of size \code{c(n, p)} with the random sample.
#' \item \code{c_ACG}: the normalizing constant.
#' }
#' @details
#' Due to the projection of the ACG, the shape matrix
#' \eqn{\boldsymbol{\Lambda}}{\Lambda} is only identified up to a constant,
#' that is, \eqn{\boldsymbol{\Lambda}}{\Lambda} and
#' \eqn{c\boldsymbol{\Lambda}}{c\Lambda} give the same ACG distribution.
#' Usually, \eqn{\boldsymbol{\Lambda}}{\Lambda} is normalized to have trace
#' equal to \eqn{p}.
#'
#' \code{c_ACG} is vectorized on \code{p}. If \eqn{p = 1}, then the ACG is the
#' uniform distribution in the set \eqn{\{-1, 1\}}{{-1, 1}}.
#' @author Eduardo García-Portugués, Davy Paindaveine, and Thomas Verdebout.
#' @references
#' Tyler, D. E. (1987). Statistical analysis for the angular central Gaussian
#' distribution on the sphere. \emph{Biometrika}, 74(3):579--589.
#' \doi{10.1093/biomet/74.3.579}
#' @examples
#' # Simulation and density evaluation for p = 2
#' Lambda <- diag(c(5, 1))
#' n <- 1e3
#' x <- r_ACG(n = n, Lambda = Lambda)
#' col <- viridisLite::viridis(n)
#' r <- runif(n, 0.95, 1.05) # Radius perturbation to improve visualization
#' plot(r * x, pch = 16, col = col[rank(d_ACG(x = x, Lambda = Lambda))])
#'
#' # Simulation and density evaluation for p = 3
#' Lambda <- rbind(c(5, 1, 0.5),
#'                 c(1, 2, 1),
#'                 c(0.5, 1, 1))
#' x <- r_ACG(n = n, Lambda = Lambda)
#' if (requireNamespace("rgl")) {
#'   rgl::plot3d(x, col = col[rank(d_ACG(x = x, Lambda = Lambda))], size = 5)
#' }
#' @name ACG

#' @rdname ACG
#' @export
d_ACG <- function(x, Lambda, log = FALSE) {

# x as matrix
if (is.null(dim(x))) {

x <- rbind(x)

}

# Dimension
p <- ncol(x)
if (p != sqrt(length(Lambda))) {

stop("x and Lambda do not have the same dimension.")

}
if (p == 1) {

# Log-density
log_dens <- d_unif_sphere(x = x, log = TRUE)

} else {

# Detect edge case x = 0 and renormalize if necessary
x <- check_unit_norm(x = x, warnings = TRUE)

# Log-density
log_dens <- c_ACG(p = p, Lambda = Lambda, log = TRUE) -
0.5 * p * log(rowSums((x %*% solve(Lambda)) * x))

}
return(switch(log + 1, exp(log_dens), log_dens))

}

#' @rdname ACG
#' @export
c_ACG <- function(p, Lambda, log = FALSE) {

# Check symmetry of Lambda
if (!isSymmetric(Lambda, tol = sqrt(.Machine\$double.eps),
check.attributes = FALSE)) {

stop("Lambda must be a symmetric matrix")

}

# Log-determinant using Cholesky decomposition (generates an error if Lambda
# is not positive definite)
log_det <- 2 * sum(log(diag(chol(Lambda))))
log_c_ACG <- - (w_p(p = p, log = TRUE) + 0.5 * log_det)
return(switch(log + 1, exp(log_c_ACG), log_c_ACG))

}

#' @rdname ACG
#' @export
r_ACG <- function(n, Lambda) {

# Dimension
p <- sqrt(length(Lambda))

# Simulation (chol() generates an error if Lambda is not positive definite)
x <- matrix(rnorm(n = n * p), nrow = n, ncol = p, byrow = TRUE) %*%
chol(Lambda)

# Projection
return(x / sqrt(rowSums(x * x)))

}


## Try the rotasym package in your browser

Any scripts or data that you put into this service are public.

rotasym documentation built on Oct. 18, 2021, 5:09 p.m.