R/uniformly_functions.R

In thames: Truncated Harmonic Mean Estimator of the Marginal Likelihood

#' Uniform sampling on/in ellipsoid
#' @description Uniform sampling on an ellipsoid or in an ellipsoid.
#'   The sampling \emph{in} an ellipsoid is available in arbitrary
#'   dimension. The sampling \emph{on} an ellipsoid is available only in
#'   dimension 2 or 3.
#'   This code is borrowed directly from the uniformly package
#'   (https://CRAN.R-project.org/package=uniformly)
#'
#' @name runif_ellipsoid
#' @rdname runif_ellipsoid
#'
#' @param n number of simulations
#' @param A symmetric positive-definite matrix defining the ellipsoid (see
#'   Details), of size 2 for \code{runif_on_ellipse} and size 2 or 3 for
#'   \code{runif_on_ellipsoid} (for size 2 these are the same functions)
#' @param r "radius" (see Details)
#'
#' @return The simulations in a matrix with \code{n} rows.
#'
#' @details The ellipsoid is the set of vectors \code{x} satisfying
#'   \code{t(x) \%*\% A \%*\% x == r^2}. For example, for an axis-aligned
#'   ellipse with horizontal radius \code{a} and vertical radius \code{b}, take
#'   \code{A=1/diag(c(a^2,b^2))} and \code{r=1}.
#'
NULL

#' @keywords internal
#' @rdname runif_ellipsoid
runif_on_ellipse <- function(n, A, r){
  stopifnot(r > 0)
  if(!is.matrix(A) || nrow(A) != 2L || ncol(A) != 2L) {
    stop(
      "The matrix `A` should be square of size 2."
    )
  }
  S <- A / (r * r)
  stopifnot(isSymmetric(S))
  e <- eigen(S, symmetric = TRUE)
  eigenvalues <- e$values
  if(any(eigenvalues <= 0)) stop("`A` is not positive.")
  ainv2 <- eigenvalues[1L]; binv2 <- eigenvalues[2L]
  a2 <- 1/ainv2; b2 <- 1/binv2
  a <- sqrt(a2); b <- sqrt(b2)
  if(isTRUE(all.equal(a, b))) {
    return(runif_on_sphere(n, 2L, a))
  }
  rotMatrix <- e$vectors
  ssims <- runif_on_sphere(n, 2L, r = 1)
  allsims <- cbind(a*ssims[, 1L], b*ssims[, 2L])
  p <- a * sqrt(binv2*sq(ssims[, 2L]) + ainv2*sq(ssims[, 1L]))
  accept <- which(runif(n) < p)
  sims <- allsims[accept, ]
  nsims <- length(accept)
  while(nsims < n) {
    newn <- n - nsims
    ssims <- runif_on_sphere(newn, 2L, r = 1)
    allsims <- cbind(a*ssims[, 1L], b*ssims[, 2L])
    p <- a * sqrt(binv2*sq(ssims[, 2L]) + ainv2*sq(ssims[, 1L]))
    accept <- which(runif(newn) < p)
    sims <- rbind(sims, allsims[accept, ])
    nsims <- nsims + length(accept)
  }
  t(rotMatrix %*% t(sims))
}

#' @keywords internal
sq <- function(x){
  x^2
}

#' @keywords internal
#' @rdname runif_ellipsoid
runif_on_ellipsoid <- function(n, A, r){
  if(is.matrix(A) && nrow(A) == 2L && ncol(A) == 2L) {
    return(runif_on_ellipse(n, A, r))
  }
  stopifnot(r > 0)
  if(!is.matrix(A) || nrow(A) != 3L || ncol(A) != 3L) {
    stop(
      "This function should be only used for sampling on the boundary of an ",
      "ellipse or a triaxial ellipsoid only."
    )
  }
  S <- A / (r * r)
  stopifnot(isSymmetric(S))
  e <- eigen(S, symmetric = TRUE)
  eigenvalues <- e$values
  if(any(eigenvalues <= 0)) stop("`A` is not positive.")
  ainv2 <- eigenvalues[1L]; binv2 <- eigenvalues[2L]; cinv2 <- eigenvalues[3L]
  a2 <- 1/ainv2; b2 <- 1/binv2; c2 <- 1/cinv2
  a <- sqrt(a2); b <- sqrt(b2); c <- sqrt(c2)
  if(isTRUE(all.equal(a, c))) {
    return(runif_on_sphere(n, 3L, a))
  }
  rotMatrix <- e$vectors
  ssims <- runif_on_sphere(n, 3L, r = 1)
  allsims <- cbind(a*ssims[, 1L], b*ssims[, 2L], c*ssims[, 3L])
  p <- a * sqrt(
    binv2*sq(ssims[, 2L]) + cinv2*sq(ssims[, 3L]) + ainv2*sq(ssims[, 1L])
  )
  accept <- which(runif(n) < p)
  sims <- allsims[accept, ]
  nsims <- length(accept)
  while(nsims < n) {
    newn <- n - nsims
    ssims <- runif_on_sphere(newn, 3L, r = 1)
    allsims <- cbind(a*ssims[, 1L], b*ssims[, 2L], c*ssims[, 3L])
    p <- a * sqrt(
      binv2*sq(ssims[, 2L]) + cinv2*sq(ssims[, 3L]) + ainv2*sq(ssims[, 1L])
    )
    accept <- which(runif(newn) < p)
    sims <- rbind(sims, allsims[accept, ])
    nsims <- nsims + length(accept)
  }
  t(rotMatrix %*% t(sims))
}

#' @keywords internal
#' @rdname runif_ellipsoid
runif_in_ellipsoid <- function(n, A, r){
  U <- chol(A)
  x <- runif_in_sphere(n, d=ncol(A), r=r)
  t(backsolve(U, t(x)))
}

#' Uniform sampling on/in sphere
#' @description Uniform sampling on a sphere or in a sphere, in arbitrary
#' dimension.
#' This code is borrowed directly from the uniformly package
#' (https://CRAN.R-project.org/package=uniformly)
#'
#' @name runif_sphere
#' @rdname runif_sphere
#'
#' @param n number of simulations
#' @param d dimension of the space
#' @param r radius of the sphere
#'
#' @return The simulations in a \code{n} times \code{d} matrix.
#' @importFrom stats runif rnorm
NULL

#' @keywords internal
#' @rdname runif_sphere
runif_on_sphere <- function(n, d, r=1){
  sims <- matrix(rnorm(n*d), nrow=n, ncol=d)
  r * sims / sqrt(apply(sims, 1L, crossprod))
}

#' @keywords internal
#' @rdname runif_sphere
runif_in_sphere <- function(n, d, r=1){
  r * runif_on_sphere(n, d, runif(n)^(1/d))
}