RTMB (version 1.8)

Simulation: Simulation

Description

An RTMB objective function can be run in 'simulation mode' where standard likelihood evaluation is replaced by corresponding random number generation. This facilitates automatic simulation under some restrictions, notably regarding the order of likelihood accumulation in the user template - see details. Simulations can be obtained directly from the model object by obj$simulate() or used indirectly via checkConsistency. In both cases a simulation is carried out from the full model i.e. both random effects and observations are simulated. The OBS function is used to mark which data items are observations - see examples.

Usage

simref(n)
# S3 method for simref
dim(x) <- value
# S3 method for simref
length(x)
# S3 method for simref
dim(x)
# S3 method for simref
is.array(x)
# S3 method for simref
is.matrix(x)
# S3 method for simref
as.array(x, ...)
# S3 method for simref
is.na(x)
# S3 method for simref
[(x, ...)
# S3 method for simref
[(x, ...) <- value
# S3 method for simref
Ops(e1, e2)
# S3 method for simref
Math(x, ...)
# S3 method for simref
t(x)
# S3 method for simref
diff(x, lag = 1L, differences = 1L, ...)
# S3 method for simref
Summary(..., na.rm = FALSE)

Value

An object with write access to store the simulation.

Arguments

n

Length

x

Object of class 'simref'

value

Replacement (numeric)

...

Extra arguments

e1

First argument

e2

Second argument

lag

As diff

differences

As diff

na.rm

Ignored

Functions

• simref(): Construct simref

• dim(simref) <- value: Equivalent of dim<-

• length(simref): Equivalent of length

• dim(simref): Equivalent of dim

• is.array(simref): Equivalent of is.array

• is.matrix(simref): Equivalent of is.matrix

• as.array(simref): Equivalent of as.array

• is.na(simref): Equivalent of is.na

• [: Equivalent of [

• `[`(simref) <- value: Equivalent of [<-

• Ops(simref): Equivalent of Ops

• Math(simref): Equivalent of Math

• t(simref): Equivalent of t

• diff(simref): Equivalent of diff

• Summary(simref): Summary operations are not invertible and will throw an error.

Details

In simulation mode all log density evaluation, involving either random effects or observations, is interpreted as probability assignment.

direct vs indirect Assignments can be 'direct' as for example

dnorm(u, log=TRUE) ## u ~ N(0, 1)

or 'indirect' as in

dnorm(2*(u+1), log=TRUE) ## u ~ N(-1, .25)

Indirect assignment works for a limited set of easily invertible functions - see methods(class="simref").

Simulation order Note that probability assignments are sequential: All information required to draw a new variable must already be simulated. It follows that, for the simulation to work, one cannot assume likelihood accumulation is commutative!

Vectorized assignment implicitly occurs elementwise from left to right. For example the assignment

dnorm(diff(u), log=TRUE)

is not valid without a prior assignment of u[1], e.g.

dnorm(u[1], log=TRUE)

Supported distributions Assignment must use supported density functions. I.e.

dpois(N, exp(u), log=TRUE)

cannot be replaced by

N * u - exp(u)

The latter will have no effect in simulation mode (the simulation will be NA).

Return value Note that when in simulation mode, the density functions all return zero. The actual simulation is written to the input argument by reference. This is very unlike standard R semantics.

Examples

## Basic example simulating response marked by OBS()
func <- function(par) {
  getAll(par, tmbdata)
  y <- OBS(y)
  ans <- -sum(dbinom(y, prob = plogis(a+b*x), size = 1, log = TRUE))
}
set.seed(101)
tmbdata <- list(x = seq(-3, 3, length = 25))
tmbdata$y <- rbinom(25, size = 1, prob = plogis(2 - 0.1*tmbdata$x))
obj <- MakeADFun(func, list(a = 0, b = 0), silent=TRUE)
with(obj, nlminb(par, fn, gr))
obj$simulate()
## Basic example simulating random effects
func <- function(p) {
  u <- p$u
  ans <- -dnorm(u[1], log=TRUE) ## u[1] ~ N(0,1)
  ans <- ans - sum(dnorm(diff(u), log=TRUE)) ## u[i]-u[i-1] ~ N(0,1)
}
obj <- MakeADFun(func, list(u=numeric(20)), random="u")
obj$simulate()
## Demonstrate how a 'simref' object works
s <- simref(4)
s2 <- 2 * s[1:2] + 1
s2[] <- 7
s ## 3 3 NA NA