R/sobol_matrices.R

Defines functions sobol_matrices scrambled_cluster scrambled_third scrambled_second scrambled_sobol

Documented in sobol_matrices

# INTERNAL FUNCTION TO CREATE THE SCRAMBLED MATRIX FOR THE
# COMPUTATION OF FIRST-ORDER SOBOL' INDICES -----------------------------------

scrambled_sobol <- function(A, B) {
  X <- rbind(A, B)
  for(i in 1:ncol(A)) {
    AB <- A
    AB[, i] <- B[, i]
    X <- rbind(X, AB)
  }
  AB <- X
  return(AB)
}

# INTERNAL FUNCTION TO CREATE THE SCRAMBLED MATRIX FOR THE
# COMPUTATION OF SECOND-ORDER SOBOL' INDICES ----------------------------------

scrambled_second <- function(A, B) {
  X <- rbind(A, B)
  parms <- utils::combn(1:ncol(A), 2, simplify = FALSE)
  for(i in parms) {
    AB <- A
    AB[, c(i[[1]], i[[2]])] <- B[, c(i[[1]], i[[2]])]
    X <- rbind(X, AB)
  }
  n <- nrow(X)
  AB <- X[(2 * nrow(A) + 1):n, ]
  return(AB)
}

# INTERNAL FUNCTION TO CREATE THE SCRAMBLED MATRIX FOR THE
# COMPUTATION OF THIRD-ORDER SOBOL' INDICES ----------------------------------

scrambled_third <- function(A, B) {
  X <- rbind(A, B)
  parms <- utils::combn(1:ncol(A), 3, simplify = FALSE)
  for(i in parms) {
    AB <- A
    AB[, c(i[[1]], i[[2]], i[[3]])] <- B[, c(i[[1]], i[[2]], i[[3]])]
    X <- rbind(X, AB)
  }
  n <- nrow(X)
  AB <- X[(2 * nrow(A) + 1):n, ]
  return(AB)
}

# INTERNAL FUNCTION TO CREATE THE SCRAMBLED MATRIX FOR THE
# COMPUTATION OF SOBOL' INDICES OF CLUSTERS OF PARAMETERS ---------------------

scrambled_cluster <- function(A, B, cluster) {
  X <- rbind(A, B)
  for(i in cluster) {
    AB <- A
    AB[, i] <- B[, i]
    X <- rbind(X, AB)
  }
  AB <- X
  return(AB)
}

# FUNCTION TO CREATE THE SOBOL' MATRIX FOR THE COMPUTATION OF
# FIRST AND TOTAL ORDER EFFECTS - AS WELL AS SECOND AND THIRD
# ORDER EFFECTS, IF DESIRED ---------------------------------------------------

#' Creation of the sample matrices
#'
#' It creates the sample matrices to compute Sobol' first and total-order indices.
#' If needed, it also creates the sample matrices required to compute second and
#' third-order indices. It uses Sobol' quasi-random number sequences.
#'
#' @param n Integer, sample size of the Sobol' matrix.
#' @param k Integer, number of model inputs.
#' @param second Logical. If \code{second = TRUE}, it creates the scrambled
#' matrix required to compute second-order indices. Default is \code{second = FALSE}.
#' @param third Logical. If \code{third = TRUE}, it creates the scrambled
#' matrix required to compute third-order indices. Default is \code{third = FALSE}.
#' @param cluster List of vectors, each vector including the column number or the
#' column name of the parameters forming the cluster. The default is \code{cluster = NULL}
#' @seealso Check the function \code{\link{sobol}} in the package \code{randtoolbox}
#' to see how the Sobol' quasi-random number sequences are constructed.
#'
#' @return A matrix.
#' @export
#' @details If \code{cluster = NULL}, the function generates an \eqn{(n, 2k)} matrix using Sobol' quasi-random
#' number sequences. The first \emph{k}-matrix is the \strong{A} matrix and the
#' remaining \emph{k}-matrix, the \strong{B} matrix. It then generates \emph{k}
#' additional matrices {(\strong{A}^j_{\strong{B}})}, \eqn{j=1,2,...,k}, where
#' the \emph{k} matrix is composed of all columns of the \strong{A} matrix
#' except the \emph{j}-th column, which is the \emph{j} column of the
#' \strong{B} matrix. This approach leds to a total number of model runs of
#' \eqn{n(k + 2)} for first and total-order indices \insertCite{Saltelli2010a}{sensobol}.
#'
#' If a list of vectors is assigned to \code{cluster}, the output is a scrambled matrix
#' with all columns from matrix \strong{A} except those of the parameters included in the cluster, which
#' come from matrix \strong{B}.
#'
#' @importFrom Rdpack reprompt
#' @references
#' \insertAllCited{}
#' @examples
#' sobol_matrices(n = 100, k = 8, second = TRUE, third = TRUE)
sobol_matrices <- function(n, k, second = FALSE, third = FALSE, cluster = NULL) {
  # Create the Sobol quasi-random number sequence
  df <- randtoolbox::sobol(n = n,
                           dim = k * 2)
  # Create A matrix
  A <- df[, 1:k]
  # Create B matrix
  B <- df[, (k + 1) : (k * 2)]
  # Create AB matrix
  AB <- scrambled_sobol(A = A, B = B)
  if(second == TRUE & third == FALSE) {
    # Compute AB matrix for second order indices
    AB.2 <- scrambled_second(A = A, B = B)
    AB <- rbind(AB, AB.2)
  }
  if(second == FALSE & third == TRUE) {
    stop("The computation of third-order Sobol' indices requires the computation of second-order indices first")
  }
  if(second == TRUE & third == TRUE) {
    AB.2 <- scrambled_second(A = A, B = B)
    AB.3 <- scrambled_third(A = A, B = B)
    AB <- rbind(AB, AB.2, AB.3)
  }
  if(is.null(cluster) == FALSE) {
    AB <- scrambled_cluster(A = A, B = B, cluster = cluster)
  }
  return(AB)
}
arnaldpuy/sensobol documentation built on April 6, 2019, 2:18 a.m.