## ##############################################################
#'
#' @title Simulate data from array.
#' @description Simulate data (slice of) an array: Simulate n
#' observations from the array x conditional on the variables in
#' margin (a vector of indices) takes values given by margin.value
#' @author Søren Højsgaard, \email{sorenh@@math.aau.dk}
#' @name array-simulate
#'
## ##############################################################
#'
#' @param x,object An array.
#' @param nsim Number of cases to simulate.
#' @param margin,value.margin Specification of slice of array to
#' simulate from.
#' @param seed Seed to be used for random number generation.
#' @param ... Additional arguments, currently not used.
#'
#' @return A matrix.
#'
#' @note The current implementation is fragile in the sense that it is
#' not checked that the input argument \code{x} is an array.
#'
#' @keywords utilities
#' @examples
#' ## 2x2 array
#' x <- parray(c("a", "b"), levels=c(2, 2), values=1:4)
#'
#' ## Simulate from entire array
#' s <- simulateArray(x, 1000)
#' xtabs(~., as.data.frame(s))
#'
#' ## Simulate from slice defined by that dimension 1 is fixed at level 2
#' s <-simulateArray(x, 6000, 1, 2)
#' xtabs(~., as.data.frame(s))
#'
#' ## 2 x 2 x 2 array
#' x <- parray(c("a", "b", "c"), levels=c(2, 2, 2), values=1:8)
#' ## Simulate from entire array
#' s <-simulateArray(x, 36000)
#' xtabs(~., as.data.frame(s))
#'
#' ## Simulate from slice defined by that dimension 3 is fixed at level 1
#' s <-simulateArray(x, 10000, 3, 1)
#' xtabs(~., as.data.frame(s))
#'
#'
#' @export simulateArray
simulateArray <- function(x, nsim=1, margin, value.margin, seed=NULL){
if(missing(margin)) {
rhs <- NULL
lhs.dim <- dim(x)
lhs.names <- names(dimnames(x))
} else {
rhs <- margin
idx <- (1:length(dim(x)))[-rhs]
lhs.dim <- (dim(x))[idx]
lhs.names <- names(dimnames(x))[-rhs]
}
##cat(sprintf("rhs=%s, lhs.dim=%s, lhs.names=%s\n",
## toString(rhs), toString(lhs.dim), toString(lhs.names)))
llhs.dim <- length(lhs.dim)
## FIXME: simulateArray uses tableSlice(); not a big issue but still
pp <- tableSlice(x, margin=rhs, level=value.margin)
##print(pp)
set.seed(seed)
samp <- sample(length(pp), size=nsim, replace=TRUE, prob=pp)
##print(samp)
cp <- cumprod(c(1, lhs.dim[-llhs.dim]))
res <- matrix(0, nrow=nsim, ncol=llhs.dim)
for(j in 1:nsim){
res[j,] <- 1 + ((samp[j] - 1) %/% cp) %% lhs.dim
}
colnames(res) <- lhs.names
res
}
#' @export
#' @rdname array-simulate
simulate.table <- function(object, nsim=1, seed=NULL, margin, value.margin, ...){
simulateArray(object, nsim=nsim, margin=margin, value.margin=value.margin, seed=seed)
}
#' @export
#' @rdname array-simulate
simulate.xtabs <- simulate.table
#' @export
#' @rdname array-simulate
simulate.array <- simulate.table
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.