R/halton.R

In philswatton/logitPlus: What the Package Does (One Line, Title Case)

#' Halton Sequences
#'
#' @param draw Number of draws for the sequence
#' @param dim Number of dimensions to draw
#' @param prime If \code{NULL}, the function will use prime numbers in ascending order. If specified, should be an atomic vector of length \code{dim} containing prime numbers to be used for each dimension.
#' @param scramble If \code{TRUE}, randomly permute the sequences
#' @param randomise If \code{TRUE}, randomise each draw
#' @param normal If \code{TRUE}, convert the sequence to a standard normal distribution
#'
#' @return
#'
#' @references
#' \insertRef{bhat2003}{logitPlus}
#'
#' \insertRef{train2003}{logitPlus}
#'
#' \insertRef{okten2012}{logitPlus}
#'
#' @export
halton <- function(draw, dim, prime=NULL, scramble=F, randomise=F, normal=F) {

  if (is.null(prime)) {prime <- logitPlus:::prime; pckgPrime <- T} else pckgPrime <- F
  if (length(prime) < dim) stop("Length of prime is smaller than number of dimensions")
  if (!pckgPrime & length(prime) > dim) warning(paste0(length(prime), " prime numbers have been included but only ", dim, " dimensions have been specified. Using the first ", dim, " prime numbers."))

  mat <- matrix(numeric(draw*dim), draw, dim)

  halton_r <- function(draw, seq, interval, t) {
    breaks <- eval(interval)
    seq <- c(seq, rowSums(expand.grid(seq, breaks)))
    if ((length(seq) + 1) >= draw) {
      return(seq[2:(draw+1)])
    } else {
      return(halton_r(draw, seq, interval, t+1))
    }
  }

  for (i in 1:dim) {
    mat[,i] <- halton_r(draw=draw, seq=0, interval=quote(1:(prime[i]-1)/(prime[i]^t)), t=1)
  }

  if (scramble) {
    for (i in 1:dim) {
      mat[,i] <- mat[sample(draw, replace=F),i]
    }
  }

  if (randomise) {
    for (i in 1:dim) {
      mu <- runif(1, 0, 1)
      mat[,i] <- mat[,i] + 1
      mat[,i] <- ifelse(mat[,1] > 1, mat[,i]-1, mat[,i])
    }
  }

  if (normal) {
    for (i in 1:dim) {
      mat[,i] <- stats::qnorm(mat[,i])
    }
  }

  return(mat)

}