R/ACG.R

In rotasym: Tests for Rotational Symmetry on the Hypersphere

#' @title Angular central Gaussian distribution
#'
#' @description Density and simulation of the Angular Central Gaussian (ACG)
#' distribution on \eqn{\mathcal{S}^{p-1}:=\{\boldsymbol{x}\in
#' \mathbb{R}^p:\|\boldsymbol{x}\|=1\}}, \eqn{p\ge 1}. The density at
#' \eqn{\boldsymbol{x} \in \mathcal{S}^{p-1}}, \eqn{p\ge 2}, is given by
#' \deqn{c^{\mathrm{ACG}}_{p,\boldsymbol{\Lambda}} (\boldsymbol{x}'
#' \boldsymbol{\Lambda}^{-1} \boldsymbol{x})^{-p/2} \quad\mathrm{with}\quad
#' c^{\mathrm{ACG}}_{p,\boldsymbol{\Lambda}}:= 1 / (\omega_p
#' |\boldsymbol{\Lambda}|^{1/2})}
#' where \eqn{\boldsymbol{\Lambda}} is the shape matrix, a \eqn{p\times p}
#' symmetric and positive definite matrix, and \eqn{\omega_p} is the surface
#' area of \eqn{\mathcal{S}^{p-1}}.
#'
#' @inheritParams unif
#' @param Lambda the shape matrix \eqn{\boldsymbol{\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}} is only identified up to a constant, that is,
#' \eqn{\boldsymbol{\Lambda}} and \eqn{c\boldsymbol{\Lambda}} give the same ACG
#' distribution. Usually, \eqn{\boldsymbol{\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\}}.
#' @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
#' dens <- d_ACG(x = x, Lambda = Lambda)
#' plot(r * x, pch = 16, col = col[rank(dens)])
#'
#' # 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)
#' dens <- d_ACG(x = x, Lambda = Lambda)
#' scatterplot3d::scatterplot3d(x, xlim = c(-1, 1), ylim = c(-1, 1),
#'                              zlim = c(-1, 1), color = col[rank(dens)],
#'                              pch = 16, xlab = "", ylab = "", zlab = "",
#'                              angle = 20)
#' @seealso \code{\link{tangent-elliptical}}, \code{\link{unif}}.
#' @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)))

}