lava (version 1.11)

sim.default: Monte Carlo simulation

Description

Applies a function repeatedly for a specified number of replications or over a list/data.frame with plot and summary methods for summarizing the Monte Carlo experiment. Can be parallelized via the future package (use the future::plan() function).

Usage

# S3 method for default
sim(
  x = NULL,
  R = 100,
  f = NULL,
  colnames = NULL,
  seed = NULL,
  args = list(),
  iter = FALSE,
  mc.cores,
  progressr.message = NULL,
  estimate.index = 1:2,
  ...
)

Arguments

x

function or 'sim' object

R

Number of replications or data.frame with parameters

f

Optional function (i.e., if x is a matrix)

colnames

Optional column names

seed

(optional) Seed passed to set.seed() before the simulation.

args

(optional) list of named arguments passed to (mc)mapply

iter

If TRUE the iteration number is passed as first argument to (mc)mapply

mc.cores

Optional number of cores. Will use parallel::mcmapply instead of future

progressr.message

Optional message for the progressr progress-bar

estimate.index

If return object inherits from estimate then only these column indices are extracted (estimate, se, lower, upper, p-val)

...

Additional arguments to future.apply::future_mapply()

Details

To parallelize the calculation use the future::plan() function (e.g., future::plan(multisession()) to distribute the calculations over the R replications on all available cores). The output is controlled via the progressr package (e.g., progressr::handlers(global=TRUE) to enable progress information).

Each replication uses its own L'Ecuyer-CMRG random number stream derived from the current state of the random number generator (or seed) as in future.apply::future_mapply() with future.seed = TRUE (the default here). The random number generator state of the calling session is restored afterwards.

See Also

summary.sim() plot.sim() sim.lvm()

Examples

m <- lvm(y~x+e)
distribution(m,~y) <- 0
distribution(m,~x) <- dist_uniform(a=-1.1,b=1.1)
transform(m,e~x) <- function(x) (1*x^4)*rnorm(length(x),sd=1)

onerun <- function(iter=NULL,...,n=2e3,b0=1,idx=2) {
    d <- sim(m,n,p=c("y~x"=b0))
    l <- lm(y~x,d)
    res <- c(coef(summary(l))[idx,1:2],
             confint(l)[idx,],
             coef(estimate(l), mat=TRUE)[idx,2:4])
    names(res) <- c("Estimate","Model.se","Model.lo","Model.hi",
                    "Sandwich.se","Sandwich.lo","Sandwich.hi")
    res
}
val <- sim(onerun,R=10,b0=1)
val

val <- sim(val,R=40,b0=1) ## append results
summary(val,estimate=c(1,1),confint=c(3,4,6,7),true=c(1,1))

summary(val,estimate=c(1,1),se=c(2,5),names=c("Model","Sandwich"))
summary(val,estimate=c(1,1),se=c(2,5),true=c(1,1),
        names=c("Model","Sandwich"),confint=TRUE)

if (interactive()) {
    plot(val,estimate=1,c(2,5),true=1,
         names=c("Model","Sandwich"),polygon=FALSE)
    plot(val,estimate=c(1,1),se=c(2,5),main=NULL,
         true=c(1,1),names=c("Model","Sandwich"),
         line.lwd=1,col=c("gray20","gray60"),
         rug=FALSE)
    plot(val,estimate=c(1,1),se=c(2,5),true=c(1,1),
         names=c("Model","Sandwich"))
}

f <- function(a=1, b=1) {
  rep(a*b, 5)
}
R <- Expand(a=1:3, b=1:3)
sim(f, R)
sim(function(a,b) f(a,b), 3, args=c(a=5,b=5))
sim(function(iter=1,a=5,b=5) iter*f(a,b), iter=TRUE, R=5)

## Returning estimate objects with extra per-iteration information
onerun2 <- function(...) {
  x <- rnorm(50)
  y <- 0.5*x + rnorm(50)
  e <- estimate(lm(y ~ x))
  c(e, converged = 1, niter = sample(5:20, 1), drop.ic = TRUE)
  # works identically with summary.estimate objects
  # c(summary(e), converged = 1, niter = sample(5:20, 1), drop.ic = TRUE)
}
val2 <- sim(onerun2, R = 10)
val2
# extra columns are stored but excluded from default summary statistics
idx <- c("(Intercept)", "x", "niter")
summary(val2, estimate=idx, true=c(0,0.5,NA))