R/regions.R

In rotations: Working with Rotation Data

#' Confidence and credible regions for the central orientation
#'
#' Find the radius of a \eqn{100(1-\alpha)}\% confidence or credible region for
#' the central orientation based on the projected mean or median. For more on
#' the currently available methods see \code{\link{prentice}},
#' \code{\link{fisheretal}}, \code{\link{chang}}, \code{\link{zhang}} and
#' \code{\link{bayesCR}}.
#'
#' @param x \eqn{n\times p}{n-by-p} matrix where each row corresponds to a
#'   random rotation in matrix (\eqn{p=9}) or quaternion (\eqn{p=4}) form.
#' @param method character string specifying which type of interval to report,
#'   "bayes", "transformation" or "direct" based theory.
#' @param type character string, "bootstrap" or "asymptotic" are available.  For
#'   Bayes regions, give the type of likelihood: "Cayley","Mises" or "Fisher."
#' @param estimator character string either "mean" or "median."  Note that not
#'   all method/type combinations are available for both estimators.
#' @param alp the alpha level desired, e.g. 0.05 or 0.10.
#' @param ... additional arguments that are method specific.
#' @return For frequentist regions only the radius of the confidence region
#'   centered at the specified estimator is returned. For Bayes regions the
#'   posterior mode and radius of the credible region centered at that mode is
#'   returned.
#' @seealso \code{\link{bayesCR}}, \code{\link{prentice}},
#'   \code{\link{fisheretal}}, \code{\link{chang}}, \code{\link{zhang}}
#' @export
#' @examples
#' Rs <- ruars(20, rvmises, kappa = 10)
#'
#' # Compare the region sizes that are currently available
#'
#' region(Rs, method = "transformation", type = "asymptotic", estimator = "mean", alp = 0.1)
#' region(Rs, method = "transformation", type = "bootstrap", estimator = "mean",
#' alp = 0.1, symm = TRUE)
#' region(Rs, method = "direct", type = "bootstrap", estimator = "mean", alp = 0.1, m = 100)
#' region(Rs, method = "direct", type = "asymptotic", estimator = "mean", alp = 0.1)
#' \donttest{
#' region(Rs, method = "Bayes", type = "Mises", estimator = "mean",
#'        S0 = mean(Rs), kappa0 = 10, tuneS = 5000, tuneK = 1, burn_in = 1000, alp = .01, m = 5000)
#' }
region <- function(x, method, type, estimator, alp = NULL, ...) {
  UseMethod("region")
}

#' @rdname region
#' @export
region.Q4 <- function(x, method, type, estimator, alp = NULL, ...) {
  # Allow for upper case arguments
  method <- tolower(method)
  type <- tolower(type)
  # Change the previous method names to the new method names
  if (method == "moment")
    method <- "direct"
  else if (method == "eigen")
    method <- "transformation"
  if (type == "theory") type <- "asymptotic"

  # Allow for abbreviations
  method <- match.arg(method, c("direct", "transformation", "bayes"))
  type <- match.arg(
    arg = type,
    choices = c("asymptotic", "bootstrap", "mises", "cayley", "fisher", "mises")
  )
  Qs <- formatQ4(x)

  if (is.null(alp)) {
    # Take a default alpha=0.1 if no level is specified
    alp <- 0.1
    warning("No alpha-level specified, 0.1 used by default.")
  }

  if (method == "transformation" & type == "asymptotic") {
    if (estimator != "mean")
      stop("The method due to Prentice is only available for the mean estimator.")
    return(prentice.Q4(x = Qs, alp = alp))
  } else if (method == "direct" & type == "bootstrap") {
    return(zhang.Q4(x = Qs, estimator = estimator, alp = alp, ...))
  } else if (method == "transformation" & type == "bootstrap") {
    if (estimator != "mean")
      stop("The method due to Fisher et al. is only available for the mean estimator.")
    return(fisheretal.Q4(x = Qs, alp = alp, ...))
  } else if (method == "direct" & type == "asymptotic") {
    return(chang.Q4(x = Qs, estimator = estimator, alp = alp))
  } else if (method == "bayes") {
    if (estimator != "mean")
      stop("Bayes confidence regions are only available for the mean estimator.")
    return(bayesCR.Q4(x = Qs, type = type, alp = alp, ...))
  } else
    stop("Please choose a correct combination of method, type and estimator. See help file.")
}

#' @rdname region
#' @export
region.SO3 <- function(x, method, type, estimator, alp = NULL, ...) {
  # Allow for upper case arguments
  method <- tolower(method)
  type <- tolower(type)

  # Change the previous method names to the new method names
  if (method == "moment")
    method <- "direct"
  else if (method == "eigen")
    method <- "transformation"
  if (type == "theory") type <- "asymptotic"

  Rs <- formatSO3(x)

  # Allow for abbreviations
  method <- match.arg(method, c("direct", "transformation", "bayes"))
  type <- match.arg(
    arg = type,
    choices = c("asymptotic", "bootstrap", "mises", "cayley", "fisher", "mises")
  )

  if (is.null(alp)) {
    # Take a default alpha=0.1 if no level is specified
    alp <- 0.1
    warning("No alpha-level specified, 0.1 used by default.")
  }

  if (method == "transformation" & type == "asymptotic") {
    if (estimator != "mean")
      stop("The method due to Prentice is only available for the mean estimator.")
    return(prentice.SO3(x = Rs, alp = alp))
  } else if (method == "direct" & type == "bootstrap") {
    return(zhang.SO3(x = Rs, estimator = estimator, alp = alp, ...))
  } else if (method == "transformation" & type == "bootstrap") {
    if (estimator != "mean")
      stop("The method due to Fisher et al. is only available for the mean estimator.")
    return(fisheretal.SO3(x = Rs, alp = alp, ...))
  } else if (method == "direct" & type == "asymptotic") {
    return(chang.SO3(x = Rs, estimator = estimator, alp = alp))
  } else if (method == "bayes") {
    if (estimator != "mean")
      stop("Bayes confidence regions are only available for the mean estimator.")
    return(bayesCR.SO3(x = Rs, type = type, alp = alp, ...))
  } else
    stop("Please choose a correct combination of method, type and estimator.  See ?region for more details.")
}

#' Transformation based asymptotic confidence region
#'
#' Find the radius of a \eqn{100(1-\alpha)}\% confidence region for the
#' projected mean based on a result from directional statistics.
#'
#' Compute the radius of a \eqn{100(1-\alpha)}\% confidence region for the
#' central orientation based on the projected mean estimator using the method
#' due to \cite{prentice1986}.  For a rotation specific version see
#' \cite{rancourt2000}. The variability in each axis is different so each axis
#' will have its own radius.
#'
#' @param x \eqn{n\times p}{n-by-p} matrix where each row corresponds to a
#'   random rotation in matrix (\eqn{p=9}) or quaternion (\eqn{p=4}) form.
#' @param alp alpha level desired, e.g. 0.05 or 0.10.
#' @return Radius of the confidence region centered at the projected mean for
#'   each of the x-, y- and z-axes.
#' @seealso \code{\link{bayesCR}}, \code{\link{fisheretal}},
#'   \code{\link{chang}}, \code{\link{zhang}}
#' @details prentice1986, rancourt2000
#' @export
#' @examples
#' Qs<-ruars(20, rcayley, kappa = 100, space = 'Q4')
#'
#' # The prentice method can be accessed from the "region" function or the "prentice" function
#' region(Qs, method = "transformation", type = "asymptotic", alp = 0.1, estimator = "mean")
#' prentice(Qs, alp = 0.1)
prentice <- function(x, alp) {
  UseMethod("prentice")
}

#' @rdname prentice
#' @export
prentice.Q4 <- function(x, alp = NULL) {
  # This takes a sample qs and returns the radius of the confidence region
  # centered at the projected mean
  if (is.null(alp)) {
    # Take a default alpha=0.1 if no level is specified
    alp <- 0.1
    warning("No alpha-level specified, 0.1 used by default.")
  }

  Qs <- x
  n <- nrow(Qs)
  Shat <- mean(Qs)
  Phat <- pMat(Shat)
  Rhat <- Qs %*% Phat
  resids <- matrix(0, n, 3)
  VarShat <- matrix(0, 3, 3)
  resids <- 2 * Rhat[, 1] * matrix(Rhat[, 2:4], n, 3)
  VarShat <- (t(resids) %*% resids) / (n - 1)
  RtR <- t(Rhat) %*% Rhat
  Ahat <- (diag(RtR[1, 1], 3, 3) - RtR[-1, -1]) / n
  Tm <- diag(n * Ahat %*% solve(VarShat) %*% Ahat)
  sqrt(stats::qchisq((1 - alp), 3) / Tm)
}

#' @rdname prentice
#' @export
prentice.SO3 <- function(x, alp = NULL) {
  Qs <- as.Q4(x)
  prentice(Qs, alp)
}

#' M-estimator theory pivotal bootstrap confidence region
#'
#' Compute the radius of a \eqn{100(1-\alpha)}\% confidence region for the
#' central orientation based on M-estimation theory.
#'
#' Compute the radius of a \eqn{100(1-\alpha)}\% confidence region for the
#' central orientation based on the projected mean estimator using the method
#' due to Zhang & Nordman (2009) (unpublished MS thesis).  By construction each
#' axis will have the same radius so the radius reported is for all three axis.
#' A normal theory version of this procedure uses the theoretical chi-square
#' limiting distribution and is given by the \code{\link{chang}} option.  This
#' method is called "direct" because it used M-estimation theory for SO(3)
#' directly instead of relying on transforming a result from directional
#' statistics as \code{\link{prentice}} and \code{\link{fisheretal}} do.
#'
#' @param x \eqn{n\times p}{n-by-p} matrix where each row corresponds to a
#'   random rotation in matrix (\eqn{p=9}) or quaternion (\eqn{p=4}) form.
#' @param estimator character string either "mean" or "median."
#' @param alp alpha level desired, e.g. 0.05 or 0.10.
#' @param m number of replicates to use to estimate the critical value.
#' @return Radius of the confidence region centered at the specified estimator.
#' @seealso \code{\link{bayesCR}}, \code{\link{prentice}},
#'   \code{\link{fisheretal}}, \code{\link{chang}}
#' @export
#' @examples
#' Rs <- ruars(20, rcayley, kappa = 100)
#'
#' # The zhang method can be accesed from the "region" function or the "zhang" function
#' # They will be different because it is a bootstrap.
#' region(Rs, method = "direct", type = "bootstrap", alp = 0.1, estimator = "mean")
#' zhang(Rs, estimator = "mean", alp = 0.1)
zhang <- function(x, estimator, alp = NULL, m = 300) {
  UseMethod("zhang")
}

#' @rdname zhang
#' @export
zhang.SO3 <- function(x,
                      estimator,
                      alp = NULL,
                      m = 300) {
  # Rs is a n-by-9 matrix where each row is an 3-by-3 rotation matrix m is the
  # number of resamples to find q_1-a alp is the level of confidence desired,
  # e.g. 0.95 or 0.90 pivot logical; should the pivotal (T) bootstrap be used or
  # nonpivotal (F)

  Rs <- formatSO3(x)

  if (estimator == "median") {
    if (is.null(alp)) {
      # Take a default alpha=0.1 if no level is specified
      alp <- 0.1
      warning("No alpha-level specified, 0.1 used by default.")
    }

    n <- nrow(Rs)
    stats <- zhangMedianC(Rs, m)
    cdtilde <- cdfunsCSO3(Rs, median(Rs))
    rad <- sqrt(as.numeric(stats::quantile(stats, 1 - alp, na.rm = TRUE)) * cdtilde[1] / (2 * n * cdtilde[2] ^ 2))
  } else if (estimator == "mean") {
    Qs <- as.Q4(Rs)
    rad <- zhang.Q4(Qs, estimator, alp, m)
  } else
    stop("Please choose an estimator mean or median.")

  min(rad, pi)
}

#' @rdname zhang
#' @export
zhang.Q4 <- function(x,
                     estimator,
                     alp = NULL,
                     m = 300) {
  if (estimator == "mean") {
    if (is.null(alp)) {
      # Take a default alpha=0.1 if no level is specified
      alp <- .1
      warning("No alpha-level specified, 0.1 used by default.")
    }

    Qs <- formatQ4(x)
    n <- nrow(Qs)
    stats <- zhangQ4(Qs, m)
    cdhat <- cdfuns(Qs, estimator)
    rad <- sqrt(as.numeric(stats::quantile(stats, 1 - alp, na.rm = TRUE)) * cdhat$c / (2 * n * cdhat$d ^ 2))
  } else if (estimator == "median") {
    Rs <- as.SO3(Qs)
    rad <- zhang.SO3(Rs, estimator, alp, m)
  } else
    stop("Please choose an estimator mean or median.")

  min(rad, pi)
}

cdfuns <- function(Qs, estimator) {
  if (estimator == 'mean') {
    Shat <- mean(Qs)
    cd <- cdfunsC(Qs, Shat)
  } else if (estimator == 'median') {
    Shat <- median(Qs)
    cd <- cdfunsCMedian(Qs, Shat)
  } else
    stop("Please choose an estimator mean or median.")

  list(c = cd[1], d = cd[2])
}

#' M-estimator asymptotic confidence region
#'
#' Compute the radius of a \eqn{100(1-\alpha)}\% confidence region for the
#' central orientation based on M-estimation theory.
#'
#' Compute the radius of a \eqn{100(1-\alpha)}\% confidence region for the
#' central orientation centered at the projected mean or median based on a
#' result due to \cite{chang2001} among others.  By construction each axis will
#' have the same radius so the radius reported is for all three axes.  This
#' method is called "direct" because it uses M-estimation theory for SO(3)
#' directly instead of relying on the transformation of a result from
#' directional statistics like \code{\link{prentice}} and
#' \code{\link{fisheretal}} do.
#'
#' @param x \eqn{n\times p}{n-by-p} matrix where each row corresponds to a
#'   random rotation in matrix (\eqn{p=9}) or quaternion (\eqn{p=4}) form.
#' @param estimator character string either "mean" or "median."
#' @param alp alpha level desired, e.g. 0.05 or 0.10.
#' @return Radius of the confidence region centered at the specified estimator.
#' @details chang2001
#' @seealso \code{\link{bayesCR}}, \code{\link{prentice}},
#'   \code{\link{fisheretal}}, \code{\link{zhang}}
#' @export
#' @examples
#' Rs <- ruars(20, rcayley, kappa = 100)
#'
#' # The chang method can be accesed from the "region" function or the "chang" function
#' region(Rs, method = "direct", type = "asymptotic", alp = 0.1, estimator = "mean")
#' chang(Rs, estimator = "mean", alp = 0.1)
chang <- function(x, estimator, alp = NULL) {
  UseMethod("chang")
}

#' @rdname chang
#' @export
chang.SO3 <- function(x, estimator, alp = NULL) {
  # Rs is a n-by-9 matrix where each row is an 3-by-3 rotation matrix
  # alp is the level of confidence desired, e.g. 0.95 or 0.90
  # pivot logical; should the pivotal (T) bootstrap be used or nonpivotal (F)

  Rs <- formatSO3(x)
  Qs <- as.Q4(Rs)
  return(chang.Q4(Qs, estimator, alp))
}

#' @rdname chang
#' @export
chang.Q4 <- function(x, estimator, alp = NULL) {
  if (is.null(alp)) {
    # Take a default alpha=0.1 if no level is specified
    alp <- .1
    warning("No alpha-level specified, 0.1 used by default.")
  }

  Qs <- formatQ4(x)
  n <- nrow(Qs)
  cdhat <- cdfuns(Qs, estimator)
  rad <- sqrt(as.numeric(stats::qchisq(1 - alp, 3)) * cdhat$c / (2 * n * cdhat$d^2))
  min(rad, pi)
}

#' Transformation based pivotal bootstrap confidence region
#'
#' Find the radius of a \eqn{100(1-\alpha)}\% confidence region for the central
#' orientation based on  transforming a result from directional statistics.
#'
#' Compute the radius of a \eqn{100(1-\alpha)}\% confidence region for the
#' central orientation based on the projected mean estimator using the method
#' for the mean polar axis as proposed in \cite{fisher1996}.  To be able to
#' reduce their method to a radius requires the additional assumption of
#' rotational symmetry, equation (10) in \cite{fisher1996}.
#'
#' @param x \eqn{n\times p}{n-by-p} matrix where each row corresponds to a
#'   random rotation in matrix (\eqn{p=9}) or quaternion (\eqn{p=4}) form.
#' @param alp alpha level desired, e.g. 0.05 or 0.10.
#' @param boot should the bootstrap or normal theory critical value be used.
#' @param m number of bootstrap replicates to use to estimate critical value.
#' @param symm logical; if TRUE (default), a symmetric region is constructed.
#' @return Radius of the confidence region centered at the projected mean.
#' @seealso \code{\link{bayesCR}}, \code{\link{prentice}}, \code{\link{chang}},
#'   \code{\link{zhang}}
#' @details fisher1996
#' @export
#' @examples
#' Qs<-ruars(20, rcayley, kappa = 100, space = 'Q4')
#'
#' # The Fisher et al. method can be accesed from the "region" function or the "fisheretal" function
#' region(Qs, method = "transformation", type = "bootstrap", alp = 0.1,
#' symm = TRUE, estimator = "mean")
#' fisheretal(Qs, alp = 0.1, boot = TRUE, symm = TRUE)
fisheretal <- function(x,
                       alp = NULL,
                       boot = TRUE,
                       m = 300,
                       symm = TRUE) {
  UseMethod("fisheretal")
}

#' @rdname fisheretal
#' @export
fisheretal.Q4 <- function(x,
                          alp = NULL,
                          boot = TRUE,
                          m = 300,
                          symm = TRUE) {
  if (is.null(alp)) {
    # Take a default alpha=0.1 if no level is specified
    alp <- .1
    warning("No alpha-level specified, 0.1 used by default.")
  }

  Qs <- formatQ4(x)

  if (boot) {
    Tstats <- fisherBootC(Qs, m, symm)
    qhat <- as.numeric(stats::quantile(Tstats, 1 - alp, na.rm = TRUE))
  } else
    qhat <- stats::qchisq(1 - alp, 3)

  rsym <- stats::optim(
    par = .05,
    fn = optimAxis,
    Qs = Qs,
    cut = qhat,
    symm = TRUE,
    method = 'Brent',
    lower = 0,
    upper = pi
  )$par

  min(rsym, pi)
}

optimAxis <- function(r, Qs, cut, symm) {
  Shat <- as.Q4(mis.axis(mean(Qs)), r)
  if (symm)
    Tm <- fisherAxisC(Qs, Shat)
  else
    Tm <- fisherAxisCSymmetric(Qs, Shat)
  (Tm - cut)^2
}

#' @rdname fisheretal
#' @export
fisheretal.SO3 <- function(x,
                           alp = NULL,
                           boot = TRUE,
                           m = 300,
                           symm = TRUE) {
  Qs <- as.Q4(x)
  fisheretal(Qs, alp, boot, m, symm)
}