Movatterモバイル変換

Type:	Package
Title:	'R' Bindings for 'TMB'
Version:	1.8
Date:	2025-10-13
Author:	Kasper Kristensen [aut, cre]
Maintainer:	Kasper Kristensen <kaskr@dtu.dk>
Description:	Native 'R' interface to 'TMB' (Template Model Builder) so models can be written entirely in 'R' rather than 'C++'. Automatic differentiation, to any order, is available for a rich subset of 'R' features, including linear algebra for dense and sparse matrices, complex arithmetic, Fast Fourier Transform, probability distributions and special functions. 'RTMB' provides easy access to model fitting and validation following the principles of Kristensen, K., Nielsen, A., Berg, C. W., Skaug, H., & Bell, B. M. (2016) <doi:10.18637/jss.v070.i05> and Thygesen, U.H., Albertsen, C.M., Berg, C.W. et al. (2017) <doi:10.1007/s10651-017-0372-4>.
License:	GPL-2 |GPL-3 [expanded from: GPL (≥ 2)]
Imports:	Rcpp (≥ 1.0.9), Matrix, methods, TMB (≥ 1.9.7), MASS
LinkingTo:	Rcpp, TMB, RcppEigen
VignetteBuilder:	knitr, rmarkdown
Suggests:	knitr, rmarkdown, igraph, tinytest, numDeriv
URL:	https://github.com/kaskr/RTMB
BugReports:	https://github.com/kaskr/RTMB/issues
NeedsCompilation:	yes
Encoding:	UTF-8
Packaged:	2025-10-13 17:30:11 UTC; kaskr
Repository:	CRAN
Date/Publication:	2025-10-14 05:10:56 UTC

RTMB: R bindings for TMB

Description

The RTMB package provides a native R interface fora subset ofTMB so you can avoid coding in C++. RTMB only affects theTMB functionMakeADFun that builds the objective function. OnceMakeADFun has been invoked, everything else isexactly the sameandmodels run as fast as if coded in C++.

Details

RTMB offers a greatly simplified interface to TMB. The TMB objective function can now be written entirely in R rather than C++ (TMB-interface). In addition, we highlight two new simplifications:

1. For most cases, simulation testing can be carried outautomatically without the need to add simulation blocks (Simulation).

2. Quantile residuals can be obtained without any essential modifications to the objective function (OSA-residuals).

The introduction vignette describes these basic features - seevignette("RTMB-introduction").

In addition to the usualMakeADFun interface, RTMB offers a lower level interface to the AD machinery (MakeTape).MakeTape replaces the functionality you would normally get in TMB using C++ functors, such as calculating derivatives inside the objective function.

The advanced vignette covers these topics - seevignette("RTMB-advanced").

Note

RTMB relies heavily on the new AD framework 'TMBad' without which this interface would not be possible.

Author(s)

Kasper Kristensen

Maintainer:kaskr@dtu.dk

Distributional assignment operator

Description

Distributional assignment operator

Usage

x %~% distr

Arguments

Details

Provides a slightly simplified syntaxinspired by, butnot compatible with, other probabilistic programming languages (e.g. BUGS/JAGS):

• x %~% distribution(...) is syntactic sugar for.nll <- .nll - sum(distribution(x,...,log=TRUE))

• The variable.nll is automatically initialized to0 and returned on exit.

Value

The updated value of the hidden variable.nll.

Note

If the shorter name~ is preferred, it can be locally overloaded using"~" <- RTMB::"%~%".

Examples

f <- function(parms) {  getAll(parms)  x %~% dnorm(mu, 1)  y %~% dpois(exp(x))}p <- list(mu=0, x=numeric(10))y <- 1:10obj <- MakeADFun(f, p, random="x")

Convert R object to AD

Description

Signify that this object should be given an AD interpretation if evaluated in an active AD context. Otherwise, keep object as is.

Usage

AD(x, force = FALSE)

Arguments

Details

AD is a generic constructor, converting plain R structures to RTMB objects if in an autodiff context. Otherwise, it does nothing (and adds virtually no computational overhead).

AD knows the following R objects:

• Numeric objects frombase, such asnumeric(),matrix(),array(), are converted to classadvector with other attributes kept intact.

• Complex objects frombase, such ascomplex(), are converted to classadcomplex.

• Sparse matrices fromMatrix, such asMatrix(),Diagonal(), are converted toadsparse.

AD provides a reliable way to avoid problems with method dispatch when mixing operand types. For instance, sub assigningx[i] <- y may be problematic whenx is numeric andy isadvector. A prior statementx <- AD(x) solves potential method dispatch issues and can therefore be used as a reliable alternative toADoverload.

Examples

## numeric object to ADAD(numeric(4), force=TRUE)## complex object to ADAD(complex(4), force=TRUE)## Convert sparse matrices (Matrix package) to AD representationF <- MakeTape(function(x) {  M <- AD(Matrix::Matrix(0,4,4))  M[1,] <- x  D <- AD(Matrix::Diagonal(4))  D@x[] <- x  M + D}, 0)F(2)

AD apply functions

Description

Thesebase apply methods have been modified to keep the AD class attribute (which would otherwise be lost).

Usage

## S4 method for signature 'advector'apply(X, MARGIN, FUN, ..., simplify = TRUE)## S4 method for signature 'ANY'sapply(X, FUN, ..., simplify = TRUE, USE.NAMES = TRUE)

Arguments

Value

Object of class"advector" with a dimension attribute.

Functions

• apply(advector): Asapply

• sapply(ANY): Assapply

Examples

F <- MakeTape(function(x) apply(matrix(x,2,2), 2, sum), numeric(4))F$jacobian(1:4)

AD complex numbers

Description

A limited set of complex number operations can be used when constructing AD tapes. The available methods are listed in this help page.

Usage

adcomplex(real, imag = rep(advector(0), length(real)))## S3 method for class 'adcomplex'Re(z)## S3 method for class 'adcomplex'Im(z)## S4 method for signature 'adcomplex'show(object)## S3 method for class 'adcomplex'dim(x)## S3 replacement method for class 'adcomplex'dim(x) <- value## S3 method for class 'adcomplex'x[...]## S3 replacement method for class 'adcomplex'x[...] <- value## S3 method for class 'adcomplex't(x)## S3 method for class 'adcomplex'length(x)## S3 method for class 'adcomplex'Conj(z)## S3 method for class 'adcomplex'Mod(z)## S3 method for class 'adcomplex'Arg(z)## S3 method for class 'adcomplex'x + y## S3 method for class 'adcomplex'x - y## S3 method for class 'adcomplex'x * y## S3 method for class 'adcomplex'x / y## S3 method for class 'adcomplex'exp(x)## S3 method for class 'adcomplex'log(x, base)## S3 method for class 'adcomplex'sqrt(x)## S4 method for signature 'adcomplex'fft(z, inverse = FALSE)## S4 method for signature 'advector'fft(z, inverse = FALSE)## S3 method for class 'adcomplex'rep(x, ...)## S3 method for class 'adcomplex'as.vector(x, mode = "any")## S3 method for class 'adcomplex'is.matrix(x)## S3 method for class 'adcomplex'as.matrix(x, ...)## S4 method for signature 'adcomplex,ANY'x %*% y## S4 method for signature 'adcomplex,ANY'solve(a, b)## S4 method for signature 'adcomplex'colSums(x)## S4 method for signature 'adcomplex'rowSums(x)## S4 method for signature 'adcomplex,ANY,ANY'diag(x)## S4 method for signature 'advector,adcomplex'Ops(e1, e2)## S4 method for signature 'adcomplex,advector'Ops(e1, e2)

Arguments

Value

Object of class"adcomplex".

Functions

• adcomplex(): Constructadcomplex vector

• Re(adcomplex): Ascomplex

• Im(adcomplex): Ascomplex

• show(adcomplex): Print method

• dim(adcomplex): Asdim

• dim(adcomplex) <- value: Asdim

• [: As[

• `[`(adcomplex) <- value: As[<-

• t(adcomplex): Ast

• length(adcomplex): Aslength

• Conj(adcomplex): Ascomplex

• Mod(adcomplex): Ascomplex

• Arg(adcomplex): Ascomplex

• +: Ascomplex

• -: Ascomplex

• *: Ascomplex

• /: Ascomplex

• exp(adcomplex): Ascomplex

• log(adcomplex): Ascomplex

• sqrt(adcomplex): Ascomplex

• fft(adcomplex): Fast Fourier Transform equivalent tofft. Notably this is themultivariate transform whenx is an array.

• fft(advector): If real input is supplied it is first converted to complex.

• rep(adcomplex): Asrep

• as.vector(adcomplex): Apply for each of real/imag

• is.matrix(adcomplex): Apply for real

• as.matrix(adcomplex): Apply for each of real/imag

• x %*% y: Complex matrix multiply

• solve(a = adcomplex, b = ANY): Complex matrix inversion and solve

• colSums(adcomplex): Apply for each of real/imag

• rowSums(adcomplex): Apply for each of real/imag

• diag(x = adcomplex, nrow = ANY, ncol = ANY): Apply for each of real/imag

• Ops(e1 = advector, e2 = adcomplex): Mixed real/complex arithmetic

• Ops(e1 = adcomplex, e2 = advector): Mixed real/complex arithmetic

Examples

## Tape using complex operationsF <- MakeTape(function(x) {  x <- as.complex(x)  y <- exp( x * ( 1 + 2i ) )  c(Re(y), Im(y))}, numeric(1))FF(1)## Complex FFT on the tapeG <- MakeTape(function(x) sum(Re(fft(x))), numeric(3))G$simplify()G$print()

AD aware numeric constructors

Description

These base constructors have been extended to keep the AD class attribute of the data argument.

Usage

## S4 method for signature 'advector,ANY,ANY'diag(x, nrow, ncol)## S4 method for signature 'advector'matrix(data = NA, nrow = 1, ncol = 1, byrow = FALSE, dimnames = NULL)## S4 method for signature 'num.'matrix(data = NA, nrow = 1, ncol = 1, byrow = FALSE, dimnames = NULL)

Arguments

Value

Object of class"advector" with a dimension attribute.

Functions

• diag(x = advector, nrow = ANY, ncol = ANY): Equivalent ofdiag

• matrix(advector): Equivalent ofmatrix

• matrix(num.): Equivalent ofmatrix

Examples

func <- function(x) {  M <- matrix(x, 2, 2)  print(class(M))  D <- diag(x)  print(class(D))  0}invisible(func(1:4))            ## 'matrix' 'array'invisible(MakeTape(func, 1:4))  ## 'advector'

AD adjoint code from R

Description

Writing custom AD adjoint derivatives from R

Usage

ADjoint(f, df, name = NULL, complex = FALSE)

Arguments

Details

Reverse mode derivatives (adjoint code) can be implemented from R using the functionADjoint. It takes as input a function of a single argumentf(x) representing the function value, and another function ofthree argumentsdf(x, y, dy) representing the adjoint derivative wrtx defined as⁠d/dx sum( f(x) * dy )⁠. Bothy anddy have the same length asf(x). The argumenty can be assumed equal tof(x) to avoid recalculation during the reverse pass. It should be assumed that all argumentsx,y,dy are vectors without any attributesexcept for dimensions, which are stored on first evaluation. The latter is convenient when implementing matrix functions (seelogdet example).Higher order derivatives automatically work provided thatdf is composed by functions thatRTMB already knows how to differentiate.

Value

A function that allows for numeric and taped evaluation.

Complex case

The argumentcomplex=TRUE specifies that the functionsf anddf are complex differentiable (holomorphic) and that argumentsx,y anddy should be assumed complex (oradcomplex). Recall that complex differentiability is a strong condition excluding many continuous functions e.g.Re,Im,Conj (see example).

Note

ADjoint may be useful when you need a special atomic function which is not yet available inRTMB, or just to experiment with reverse mode derivatives.However, the approach may cause asignificant overhead compared to nativeRTMB derivatives. In addition, the approach isnot thread safe, i.e. calling R functions cannot be done in parallel using OpenMP.

Examples

############################################################################## Lambert W-function defined by W(y*exp(y))=yW <- function(x) {  logx <- log(x)  y <- pmax(logx, 0)  while (any(abs(logx - log(y) - y) > 1e-9, na.rm = TRUE)) {      y <- y - (y - exp(logx - y)) / (1 + y)  }  y}## DerivativesdW <- function(x, y, dy) {   dy / (exp(y) * (1. + y))}## Define new derivative symbolLamW <- ADjoint(W, dW)## Test derivatives(F <- MakeTape(function(x)sum(LamW(x)), numeric(3)))F(1:3)F$print()                ## Note the 'name'F$jacobian(1:3)          ## gradientF$jacfun()$jacobian(1:3) ## hessian############################################################################## Log determinantlogdet <- ADjoint(   function(x) determinant(x, log=TRUE)$modulus,   function(x, y, dy) t(solve(x)) * dy,   name = "logdet")(F <- MakeTape(logdet, diag(2)))## Test derivatives## Compare with numDeriv::hessian(F, matrix(1:4,2))F$jacfun()$jacobian(matrix(1:4,2)) ## Hessian############################################################################## Holomorphic extension of 'solve'matinv <- ADjoint(   solve,   function(x,y,dy) -t(y) %*% dy %*% t(y),   complex=TRUE)(F <- MakeTape(function(x) Im(matinv(x+AD(1i))), diag(2)))## Test derivatives## Compare with numDeriv::jacobian(F, matrix(1:4,2))F$jacobian(matrix(1:4,2))

AD matrix methods (sparse and dense)

Description

Matrices (base package) and sparse matrices (Matrix package) can be used inside theRTMB objective function as part of the calculations. Behind the scenes these R objects are converted to AD representations when needed. AD objects have a temporary lifetime, so you probably won't see them / need to know them. The only important thing is whichmethods work for the objects.

Usage

## S3 method for class 'advector'chol(x, ...)## S3 method for class 'advector'determinant(x, logarithm = TRUE, ...)## S4 method for signature 'adcomplex'eigen(x, symmetric, only.values = FALSE, EISPACK = FALSE)## S4 method for signature 'advector'eigen(x, symmetric, only.values = FALSE, EISPACK = FALSE)## S4 method for signature 'advector'svd(x, nu, nv, LINPACK = FALSE)## S3 method for class 'adsparse't(x)## S3 method for class 'adsparse'x[...]## S3 replacement method for class 'adsparse'x[...] <- value## S3 method for class 'adsparse'as.matrix(x, ...)## S4 method for signature 'adsparse,missing,missing'diag(x)## S4 method for signature 'adsparse'band(x, k1, k2)## S4 method for signature 'adsparse'tril(x, k)## S4 method for signature 'adsparse'triu(x, k)## S4 method for signature 'advector'expm(x)## S4 method for signature 'adsparse'expm(x)## S4 method for signature 'adsparse'dim(x)## S4 method for signature 'anysparse,ad'x %*% y## S4 method for signature 'ad,anysparse'x %*% y## S4 method for signature 'adsparse,adsparse'x %*% y## S4 method for signature 'ad,ad'x %*% y## S4 method for signature 'ad,ad.'tcrossprod(x, y)## S4 method for signature 'ad,ad.'crossprod(x, y)## S4 method for signature 'advector'cov2cor(V)## S4 method for signature 'ad,ad.'solve(a, b)## S4 method for signature 'num,num.'solve(a, b)## S4 method for signature 'anysparse,ad.'solve(a, b)## S4 method for signature 'advector'colSums(x, na.rm, dims)## S4 method for signature 'advector'rowSums(x, na.rm, dims)## S4 method for signature 'adsparse'colSums(x, na.rm, dims)## S4 method for signature 'adsparse'rowSums(x, na.rm, dims)## S3 method for class 'advector'cbind(...)## S3 method for class 'advector'rbind(...)## S4 method for signature 'adsparse'Math(x)

Arguments

Value

List (vectors/values) withadcomplex components.

List (vectors/values) withadvector components in symmetric case andadcomplex components otherwise.

Object of classadvector with a dimension attribute for dense matrix operations; Object of classadsparse for sparse matrix operations.

Functions

• chol(advector): AD matrix cholesky

• determinant(advector): AD log determinant

• eigen(adcomplex): General AD eigen decomposition for complex matrices. Note that argumentsymmetric isnot auto-detected somust be specified.

• eigen(advector): AD eigen decomposition for real matrices. The non-symmetric case is redirected to theadcomplex method. Note that argumentsymmetric isnot auto-detected somust be specified.

• svd(advector): AD svd decomposition for real matrices.

• t(adsparse): AD sparse matrix transpose. Re-directs tot,CsparseMatrix-method.

• [: AD sparse matrix subsetting. Re-directs to[-methods.

• `[`(adsparse) <- value: AD sparse matrix subset assignment. Re-directs to[<–methods.

• as.matrix(adsparse): Convert AD sparse to dense matrix.

• diag(x = adsparse, nrow = missing, ncol = missing): AD sparse matrix diagonal extract. Re-directs todiag,CsparseMatrix-method.

• band(adsparse): AD sparse matrix band extract. Re-directs toband,CsparseMatrix-method.

• tril(adsparse): AD sparse matrix lower triangle extract. Re-directs totril,CsparseMatrix-method.

• triu(adsparse): AD sparse matrix upper triangle extract. Re-directs totriu,CsparseMatrix-method.

• expm(advector): AD matrix exponential

• expm(adsparse): AD matrix exponential

• dim(adsparse): AD sparse matrix dimension

• x %*% y: AD matrix multiply

• x %*% y: AD matrix multiply

• x %*% y: AD matrix multiply

• x %*% y: AD matrix multiply

• tcrossprod(x = ad, y = ad.): AD matrix multiply

• crossprod(x = ad, y = ad.): AD matrix multiply

• cov2cor(advector): AD matrix cov2cor

• solve(a = ad, b = ad.): AD matrix inversion and solve

• solve(a = num, b = num.): AD matrix inversion and solve

• solve(a = anysparse, b = ad.): Sparse AD matrix solve

• colSums(advector): AD matrix (or array) colsums

• rowSums(advector): AD matrix (or array) rowsums

• colSums(adsparse): AD sparse matrix colsums

• rowSums(adsparse): AD sparse matrix rowsums

• cbind(advector): AD matrix column bind

• rbind(advector): AD matrix row bind

• Math(adsparse): AD sparse matrix 'Math group' works for functions that preserve sparsity.

Examples

F <- MakeTape(function(x) matrix(1:9,3,3) %*% x, numeric(3))F$jacobian(1:3)F <- MakeTape(function(x) Matrix::expm(matrix(x,2,2)), numeric(4))F$jacobian(1:4)F <- MakeTape(det, diag(2)) ## Indirectly available via 'determinant'F$jacobian(matrix(1:4,2))

Enable extra RTMB convenience methods

Description

Enable extra RTMB convenience methods

Usage

ADoverload(x = c("[<-", "c", "diag<-"))

Arguments

Details

Work around limitations in R's method dispatch system by overloading some selected primitives, currently:

• Inplace replacement, so you can dox[i] <- y whenx is numeric andy is AD.

• Mixed combine, so you can do e.g.c(x, y) whenx numeric andy is AD.

• Diagonal assignment, so you can dodiag(x) <- y whenx is a numeric matrix andy is AD.

In all cases, the result should be AD.The methods are automaticallytemporarily attached to the search path (search()) when enteringMakeTape orMakeADFun.Alternatively, methods can be overloaded locally inside functions using e.g."[<-" <- ADoverload("[<-"). This is only needed when using RTMB from a package.

Value

Function representing the overload.

Examples

MakeTape(function(x) {print(search()); x}, numeric(0))MakeTape(function(x) c(1,x), 1:3)MakeTape(function(x) {y <- 1:3; y[2] <- x; y}, 1)MakeTape(function(x) {y <- matrix(0,3,3); diag(y) <- x; y}, 1:3)

AD sparse matrix class

Description

Sparse matrices inRTMB are essentiallydgCMatrix with anadvector x-slot.

Slots

x

Non-zeros

i

row indices (zero based)

p

col pointers (zero based)

Dim

Dimension

The AD vector and its methods

Description

Anadvector is a class used behind the scenes to replacenormal R numeric objects during automatic differentiation. Anadvector has a temporary lifetime and therefore you do notsee /need to know it as a normal user.

Usage

advector(x)## S3 method for class 'advector'Ops(e1, e2)## S3 method for class 'advector'Math(x, ...)## S3 method for class 'advector'as.vector(x, mode = "any")## S3 method for class 'advector'as.complex(x, ...)## S3 method for class 'advector'aperm(a, perm, ...)## S3 method for class 'advector'c(...)## S3 method for class 'advector'x[...]## S3 replacement method for class 'advector'x[...] <- value## S3 method for class 'advector'x[[...]]## S3 method for class 'advector'rep(x, ...)## S3 method for class 'advector'is.nan(x)## S3 method for class 'advector'is.finite(x)## S3 method for class 'advector'is.infinite(x)## S3 method for class 'advector'is.na(x)## S3 method for class 'advector'sum(x, ..., na.rm = FALSE)## S3 method for class 'advector'mean(x, ...)## S3 method for class 'advector'prod(x, ..., na.rm = FALSE)## S3 method for class 'advector'min(..., na.rm = FALSE)## S3 method for class 'advector'max(..., na.rm = FALSE)## S3 method for class 'advector'is.numeric(x)## S3 method for class 'advector'as.double(x, ...)## S3 method for class 'advector'Complex(z)## S3 method for class 'advector'Summary(..., na.rm = FALSE)## S3 method for class 'advector'diff(x, lag = 1L, differences = 1L, ...)## S3 method for class 'advector'print(x, ...)## S4 method for signature 'num,ad,ad'ifelse(test, yes, no)## S4 method for signature 'num,num,num'ifelse(test, yes, no)## S4 method for signature 'advector'sort(x)## S4 method for signature 'advector,advector,ANY,ANY'x[i]## S4 method for signature 'advector,advector'findInterval(x, vec)

Arguments

Details

An AD vector (class='advector') is an atomic R vector of 'codes'that are internally interpretable as 'AD scalars'. A substantialpart of R's existing S3 matrix and array functionality can bere-used for AD vectors.

Value

Object of class"advector".

Functions

• advector(): Construct a new advector

• Ops(advector): Binary operations

• Math(advector): Unary operations

• as.vector(advector): Makesarray(x) work.

• as.complex(advector): Convert toADcomplex. Note that dimensions are dropped for consistency with base R.

• aperm(advector): Equivalent ofaperm

• c(advector): Equivalent ofc. However note the limitation for mixed types: Ifx is an AD type,c(x,1) works whilec(1,x) does not!

• [: Equivalent of[

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

• [[: Equivalent of[[

• rep(advector): Equivalent ofrep. Makesouter(x,x,...) work.

• is.nan(advector): Equivalent ofis.nan. Check NaN status of aconstantadvector expression. If not constant throw an error.

• is.finite(advector): Equivalent ofis.finite. Check finite status of aconstantadvector expression. If not constant throw an error.

• is.infinite(advector): Equivalent ofis.infinite. Check infinity status of aconstantadvector expression. If not constant throw an error.

• is.na(advector): Equivalent ofis.na. Check NA status of anadvector. NAs can only occur directly (as constants) or indirectly as the result of an operation with NA operands. For a tape built with non-NA parameters the NA status of any expression is constant and can therefore safely be used as part of the calculations. (assuming correct propagation of NAs via C-level arithmetic).

• sum(advector): Equivalent ofsum.na.rm=TRUE is allowed, but note that this feature assumes correct propagation of NAs via C-level arithmetic.

• mean(advector): Equivalent ofmean except no arguments beyondx are supported.

• prod(advector): Equivalent ofprod.

• min(advector): Equivalent ofmin.

• max(advector): Equivalent ofmax.

• is.numeric(advector): Makescov2cor() work. FIXME: Any unwanted side-effects with this?

• as.double(advector): Makesas.numeric() work.

• Complex(advector):Complex operations are redirected toadcomplex.

• Summary(advector): UnimplementedSummary operations (currentlyallanyrange) will throw an error.

• diff(advector): Equivalent ofdiff

• print(advector): Print method

• ifelse(test = num, yes = ad, no = ad): Equivalent ofifelse

• ifelse(test = num, yes = num, no = num): Default method

• sort(advector): Taped sorting of an AD vector

• x[i: Taped subsetting of an AD vector

• findInterval(x = advector, vec = advector): Taped interval finding of an AD vectorMakeTape(function(x) findInterval(x, AD(0:10)), 1:3)

Examples

x <- advector(1:9)a <- array(x, c(3,3))  ## as an arrayouter(x, x, "+") ## Implicit via 'rep'rev(x)           ## Implicit via '['MakeTape(sort, numeric(3))MakeTape(function(x) AD(rivers)[x], 1:3)

Distributions and special functions for which AD is implemented

Description

The functions listed in this help page are all applicable for AD types.Method dispatching follows a simple rule:If at least one argument is an AD type then a special ADimplementation is selected. In all other cases a defaultimplementation is used (typically that of thestatspackage).Argument recycling follows the R standard (although wihout any warnings).

Usage

## S4 method for signature 'ad,ad'besselK(x, nu, expon.scaled = FALSE)## S4 method for signature 'num,num'besselK(x, nu, expon.scaled = FALSE)## S4 method for signature 'ad,ad'besselI(x, nu, expon.scaled = FALSE)## S4 method for signature 'num,num'besselI(x, nu, expon.scaled = FALSE)## S4 method for signature 'ad,ad'besselJ(x, nu)## S4 method for signature 'num,num'besselJ(x, nu)## S4 method for signature 'ad,ad'besselY(x, nu)## S4 method for signature 'num,num'besselY(x, nu)## S4 method for signature 'ad,ad.,logical.'dexp(x, rate = 1, log = FALSE)## S4 method for signature 'num,num.,logical.'dexp(x, rate = 1, log = FALSE)## S4 method for signature 'osa,ANY,ANY'dexp(x, rate = 1, log = FALSE)## S4 method for signature 'simref,ANY,ANY'dexp(x, rate = 1, log = FALSE)## S4 method for signature 'ad,ad,ad.,logical.'dweibull(x, shape, scale = 1, log = FALSE)## S4 method for signature 'num,num,num.,logical.'dweibull(x, shape, scale = 1, log = FALSE)## S4 method for signature 'osa,ANY,ANY,ANY'dweibull(x, shape, scale = 1, log = FALSE)## S4 method for signature 'simref,ANY,ANY,ANY'dweibull(x, shape, scale = 1, log = FALSE)## S4 method for signature 'ad,ad,ad,logical.'dbinom(x, size, prob, log = FALSE)## S4 method for signature 'num,num,num,logical.'dbinom(x, size, prob, log = FALSE)## S4 method for signature 'osa,ANY,ANY,ANY'dbinom(x, size, prob, log = FALSE)## S4 method for signature 'simref,ANY,ANY,ANY'dbinom(x, size, prob, log = FALSE)## S4 method for signature 'ad,ad,ad,missing,logical.'dbeta(x, shape1, shape2, log)## S4 method for signature 'num,num,num,missing,logical.'dbeta(x, shape1, shape2, log)## S4 method for signature 'osa,ANY,ANY,ANY,ANY'dbeta(x, shape1, shape2, log)## S4 method for signature 'simref,ANY,ANY,ANY,ANY'dbeta(x, shape1, shape2, log)## S4 method for signature 'ad,ad,ad,missing,logical.'df(x, df1, df2, log)## S4 method for signature 'num,num,num,missing,logical.'df(x, df1, df2, log)## S4 method for signature 'osa,ANY,ANY,ANY,ANY'df(x, df1, df2, log)## S4 method for signature 'simref,ANY,ANY,ANY,ANY'df(x, df1, df2, log)## S4 method for signature 'ad,ad.,ad.,logical.'dlogis(x, location = 0, scale = 1, log = FALSE)## S4 method for signature 'num,num.,num.,logical.'dlogis(x, location = 0, scale = 1, log = FALSE)## S4 method for signature 'osa,ANY,ANY,ANY'dlogis(x, location = 0, scale = 1, log = FALSE)## S4 method for signature 'simref,ANY,ANY,ANY'dlogis(x, location = 0, scale = 1, log = FALSE)## S4 method for signature 'ad,ad,missing,logical.'dt(x, df, log)## S4 method for signature 'num,num,missing,logical.'dt(x, df, log)## S4 method for signature 'osa,ANY,ANY,ANY'dt(x, df, log)## S4 method for signature 'simref,ANY,ANY,ANY'dt(x, df, log)## S4 method for signature 'ad,ad,ad,missing,logical.'dnbinom(x, size, prob, log)## S4 method for signature 'num,num,num,missing,logical.'dnbinom(x, size, prob, log)## S4 method for signature 'osa,ANY,ANY,ANY,ANY'dnbinom(x, size, prob, log)## S4 method for signature 'simref,ANY,ANY,ANY,ANY'dnbinom(x, size, prob, log)## S4 method for signature 'ad,ad,logical.'dpois(x, lambda, log = FALSE)## S4 method for signature 'num,num,logical.'dpois(x, lambda, log = FALSE)## S4 method for signature 'osa,ANY,ANY'dpois(x, lambda, log = FALSE)## S4 method for signature 'simref,ANY,ANY'dpois(x, lambda, log = FALSE)## S4 method for signature 'ad,ad,missing,ad.,logical.'dgamma(x, shape, scale, log)## S4 method for signature 'num,num,missing,num.,logical.'dgamma(x, shape, scale, log)## S4 method for signature 'osa,ANY,ANY,ANY,ANY'dgamma(x, shape, scale, log)## S4 method for signature 'simref,ANY,ANY,ANY,ANY'dgamma(x, shape, scale, log)## S4 method for signature 'ad,ad.,ad.,missing,missing'pnorm(q, mean, sd)## S4 method for signature 'num,num.,num.,missing,missing'pnorm(q, mean, sd)## S4 method for signature 'ad,ad,missing,ad.,missing,missing'pgamma(q, shape, scale)## S4 method for signature 'num,num,missing,num.,missing,missing'pgamma(q, shape, scale)## S4 method for signature 'ad,ad,missing,missing'ppois(q, lambda)## S4 method for signature 'num,num,missing,missing'ppois(q, lambda)## S4 method for signature 'ad,ad.,missing,missing'pexp(q, rate)## S4 method for signature 'num,num.,missing,missing'pexp(q, rate)## S4 method for signature 'ad,ad,ad.,missing,missing'pweibull(q, shape, scale)## S4 method for signature 'num,num,num.,missing,missing'pweibull(q, shape, scale)## S4 method for signature 'ad,ad,ad,missing,missing,missing'pbeta(q, shape1, shape2)## S4 method for signature 'num,num,num,missing,missing,missing'pbeta(q, shape1, shape2)## S4 method for signature 'ad,ad.,ad.,missing,missing'qnorm(p, mean, sd)## S4 method for signature 'num,num.,num.,missing,missing'qnorm(p, mean, sd)## S4 method for signature 'ad,ad,missing,ad.,missing,missing'qgamma(p, shape, scale)## S4 method for signature 'num,num,missing,num.,missing,missing'qgamma(p, shape, scale)## S4 method for signature 'ad,ad.,missing,missing'qexp(p, rate)## S4 method for signature 'num,num.,missing,missing'qexp(p, rate)## S4 method for signature 'ad,ad,ad.,missing,missing'qweibull(p, shape, scale)## S4 method for signature 'num,num,num.,missing,missing'qweibull(p, shape, scale)## S4 method for signature 'ad,ad,ad,missing,missing,missing'qbeta(p, shape1, shape2)## S4 method for signature 'num,num,num,missing,missing,missing'qbeta(p, shape1, shape2)## S4 method for signature 'ad,ad'lbeta(a, b)## S4 method for signature 'num,num'lbeta(a, b)dbinom_robust(x, size, logit_p, log = FALSE)dsn(x, alpha, log = FALSE)dSHASHo(x, mu, sigma, nu, tau, log = FALSE)dtweedie(x, mu, phi, p, log = FALSE)dnbinom_robust(x, log_mu, log_var_minus_mu, log = FALSE)dnbinom2(x, mu, var, log = FALSE)dlgamma(x, shape, scale, log = FALSE)logspace_add(logx, logy)logspace_sub(logx, logy)## S4 method for signature 'ad,ad.,ad.,logical.'dnorm(x, mean = 0, sd = 1, log = FALSE)## S4 method for signature 'num,num.,num.,logical.'dnorm(x, mean = 0, sd = 1, log = FALSE)## S4 method for signature 'osa,ANY,ANY,ANY'dnorm(x, mean = 0, sd = 1, log = FALSE)## S4 method for signature 'simref,ANY,ANY,ANY'dnorm(x, mean = 0, sd = 1, log = FALSE)## S4 method for signature 'ANY,ANY,ANY,ANY'dlnorm(x, meanlog = 0, sdlog = 1, log = FALSE)## S4 method for signature 'osa,ANY,ANY,ANY'dlnorm(x, meanlog = 0, sdlog = 1, log = FALSE)## S4 method for signature 'num,num.,num.,logical.'dlnorm(x, meanlog = 0, sdlog = 1, log = FALSE)## S4 method for signature 'advector,missing,missing,missing,missing'plogis(q)## S4 method for signature 'advector,missing,missing,missing,missing'qlogis(p)dcompois(x, mode, nu, log = FALSE)dcompois2(x, mean, nu, log = FALSE)## S4 method for signature 'ad,ad,ad,missing,missing'pbinom(q, size, prob)## S4 method for signature 'num,num,num,missing,missing'pbinom(q, size, prob)## S4 method for signature 'ad,ad.,ad,logical.'dmultinom(x, size = NULL, prob, log = FALSE)## S4 method for signature 'num,num.,num,logical.'dmultinom(x, size = NULL, prob, log = FALSE)## S4 method for signature 'osa,ANY,ANY,ANY'dmultinom(x, size = NULL, prob, log = FALSE)## S4 method for signature 'simref,ANY,ANY,ANY'dmultinom(x, size = NULL, prob, log = FALSE)## S4 method for signature 'ANY,ANY,ANY,ANY'dmultinom(x, size = NULL, prob, log = FALSE)## S4 method for signature 'ad,ad.,ad.,logical.'dcauchy(x, location = 0, scale = 1, log = FALSE)## S4 method for signature 'num,num.,num.,logical.'dcauchy(x, location = 0, scale = 1, log = FALSE)## S4 method for signature 'osa,ANY,ANY,ANY'dcauchy(x, location = 0, scale = 1, log = FALSE)## S4 method for signature 'simref,ANY,ANY,ANY'dcauchy(x, location = 0, scale = 1, log = FALSE)## S4 method for signature 'ad,ad,ad,missing,logical.'dgamma(x, shape, rate, log)## S4 method for signature 'ad,ad,ad,missing,missing,missing'pnbinom(q, size, prob)

Arguments

Details

Specific documentation of the functions and arguments should be looked up elsewhere:

• All S4 methods behave as the corresponding functions in thestats package. However, some arguments may not beimplemented in the AD case (e.g.lower-tail).

• Other funtions behave as the corresponding TMB versions forwhich documentation should be looked up online.

Value

In autodiff contexts an object of class"advector" is returned; otherwise a standard numeric vector.

Functions

• besselK(x = ad, nu = ad): AD implementation ofbesselK

• besselK(x = num, nu = num): Default method

• besselI(x = ad, nu = ad): AD implementation ofbesselI

• besselI(x = num, nu = num): Default method

• besselJ(x = ad, nu = ad): AD implementation ofbesselJ

• besselJ(x = num, nu = num): Default method

• besselY(x = ad, nu = ad): AD implementation ofbesselY

• besselY(x = num, nu = num): Default method

• dexp(x = ad, rate = ad., log = logical.): AD implementation ofdexp

• dexp(x = num, rate = num., log = logical.): Default method

• dexp(x = osa, rate = ANY, log = ANY): OSA implementation

• dexp(x = simref, rate = ANY, log = ANY): Simulation implementation. Modifiesx and returns zero.

• dweibull(x = ad, shape = ad, scale = ad., log = logical.): AD implementation ofdweibull

• dweibull(x = num, shape = num, scale = num., log = logical.): Default method

• dweibull(x = osa, shape = ANY, scale = ANY, log = ANY): OSA implementation

• dweibull(x = simref, shape = ANY, scale = ANY, log = ANY): Simulation implementation. Modifiesx and returns zero.

• dbinom(x = ad, size = ad, prob = ad, log = logical.): AD implementation ofdbinom

• dbinom(x = num, size = num, prob = num, log = logical.): Default method

• dbinom(x = osa, size = ANY, prob = ANY, log = ANY): OSA implementation

• dbinom(x = simref, size = ANY, prob = ANY, log = ANY): Simulation implementation. Modifiesx and returns zero.

• dbeta(x = ad, shape1 = ad, shape2 = ad, ncp = missing, log = logical.): AD implementation ofdbeta

• dbeta(x = num, shape1 = num, shape2 = num, ncp = missing, log = logical.): Default method

• dbeta(x = osa, shape1 = ANY, shape2 = ANY, ncp = ANY, log = ANY): OSA implementation

• dbeta(x = simref, shape1 = ANY, shape2 = ANY, ncp = ANY, log = ANY): Simulation implementation. Modifiesx and returns zero.

• df(x = ad, df1 = ad, df2 = ad, ncp = missing, log = logical.): AD implementation ofdf

• df(x = num, df1 = num, df2 = num, ncp = missing, log = logical.): Default method

• df(x = osa, df1 = ANY, df2 = ANY, ncp = ANY, log = ANY): OSA implementation

• df(x = simref, df1 = ANY, df2 = ANY, ncp = ANY, log = ANY): Simulation implementation. Modifiesx and returns zero.

• dlogis(x = ad, location = ad., scale = ad., log = logical.): AD implementation ofdlogis

• dlogis(x = num, location = num., scale = num., log = logical.): Default method

• dlogis(x = osa, location = ANY, scale = ANY, log = ANY): OSA implementation

• dlogis(x = simref, location = ANY, scale = ANY, log = ANY): Simulation implementation. Modifiesx and returns zero.

• dt(x = ad, df = ad, ncp = missing, log = logical.): AD implementation ofdt

• dt(x = num, df = num, ncp = missing, log = logical.): Default method

• dt(x = osa, df = ANY, ncp = ANY, log = ANY): OSA implementation

• dt(x = simref, df = ANY, ncp = ANY, log = ANY): Simulation implementation. Modifiesx and returns zero.

• dnbinom(x = ad, size = ad, prob = ad, mu = missing, log = logical.): AD implementation ofdnbinom

• dnbinom(x = num, size = num, prob = num, mu = missing, log = logical.): Default method

• dnbinom(x = osa, size = ANY, prob = ANY, mu = ANY, log = ANY): OSA implementation

• dnbinom(x = simref, size = ANY, prob = ANY, mu = ANY, log = ANY): Simulation implementation. Modifiesx and returns zero.

• dpois(x = ad, lambda = ad, log = logical.): AD implementation ofdpois

• dpois(x = num, lambda = num, log = logical.): Default method

• dpois(x = osa, lambda = ANY, log = ANY): OSA implementation

• dpois(x = simref, lambda = ANY, log = ANY): Simulation implementation. Modifiesx and returns zero.

• dgamma(x = ad, shape = ad, rate = missing, scale = ad., log = logical.): AD implementation ofdgamma

• dgamma(x = num, shape = num, rate = missing, scale = num., log = logical.): Default method

• dgamma(x = osa, shape = ANY, rate = ANY, scale = ANY, log = ANY): OSA implementation

• dgamma(x = simref, shape = ANY, rate = ANY, scale = ANY, log = ANY): Simulation implementation. Modifiesx and returns zero.

• pnorm(q = ad, mean = ad., sd = ad., lower.tail = missing, log.p = missing): AD implementation ofpnorm

• pnorm(q = num, mean = num., sd = num., lower.tail = missing, log.p = missing): Default method

• pgamma( q = ad, shape = ad, rate = missing, scale = ad., lower.tail = missing, log.p = missing): AD implementation ofpgamma

• pgamma( q = num, shape = num, rate = missing, scale = num., lower.tail = missing, log.p = missing): Default method

• ppois(q = ad, lambda = ad, lower.tail = missing, log.p = missing): AD implementation ofppois

• ppois(q = num, lambda = num, lower.tail = missing, log.p = missing): Default method

• pexp(q = ad, rate = ad., lower.tail = missing, log.p = missing): AD implementation ofpexp

• pexp(q = num, rate = num., lower.tail = missing, log.p = missing): Default method

• pweibull( q = ad, shape = ad, scale = ad., lower.tail = missing, log.p = missing): AD implementation ofpweibull

• pweibull( q = num, shape = num, scale = num., lower.tail = missing, log.p = missing): Default method

• pbeta( q = ad, shape1 = ad, shape2 = ad, ncp = missing, lower.tail = missing, log.p = missing): AD implementation ofpbeta

• pbeta( q = num, shape1 = num, shape2 = num, ncp = missing, lower.tail = missing, log.p = missing): Default method

• qnorm(p = ad, mean = ad., sd = ad., lower.tail = missing, log.p = missing): AD implementation ofqnorm

• qnorm(p = num, mean = num., sd = num., lower.tail = missing, log.p = missing): Default method

• qgamma( p = ad, shape = ad, rate = missing, scale = ad., lower.tail = missing, log.p = missing): AD implementation ofqgamma

• qgamma( p = num, shape = num, rate = missing, scale = num., lower.tail = missing, log.p = missing): Default method

• qexp(p = ad, rate = ad., lower.tail = missing, log.p = missing): AD implementation ofqexp

• qexp(p = num, rate = num., lower.tail = missing, log.p = missing): Default method

• qweibull( p = ad, shape = ad, scale = ad., lower.tail = missing, log.p = missing): AD implementation ofqweibull

• qweibull( p = num, shape = num, scale = num., lower.tail = missing, log.p = missing): Default method

• qbeta( p = ad, shape1 = ad, shape2 = ad, ncp = missing, lower.tail = missing, log.p = missing): AD implementation ofqbeta

• qbeta( p = num, shape1 = num, shape2 = num, ncp = missing, lower.tail = missing, log.p = missing): Default method

• lbeta(a = ad, b = ad): AD implementation oflbeta

• lbeta(a = num, b = num): Default method

• dbinom_robust(): AD implementation

• dsn(): AD implementation

• dSHASHo(): AD implementation

• dtweedie(): AD implementation

• dnbinom_robust(): AD implementation

• dnbinom2(): AD implementation

• dlgamma(): AD implementation

• logspace_add(): AD implementation

• logspace_sub(): AD implementation

• dnorm(x = ad, mean = ad., sd = ad., log = logical.): AD implementation ofdnorm

• dnorm(x = num, mean = num., sd = num., log = logical.): Default method

• dnorm(x = osa, mean = ANY, sd = ANY, log = ANY): OSA implementation

• dnorm(x = simref, mean = ANY, sd = ANY, log = ANY): Simulation implementation. Modifiesx and returns zero.

• dlnorm(x = ANY, meanlog = ANY, sdlog = ANY, log = ANY): AD implementation ofdlnorm.

• dlnorm(x = osa, meanlog = ANY, sdlog = ANY, log = ANY): OSA implementation.

• dlnorm(x = num, meanlog = num., sdlog = num., log = logical.): Default method.

• plogis( q = advector, location = missing, scale = missing, lower.tail = missing, log.p = missing): Minimal AD implementation ofplogis

• qlogis( p = advector, location = missing, scale = missing, lower.tail = missing, log.p = missing): Minimal AD implementation ofqlogis

• dcompois(): Conway-Maxwell-Poisson. Calculate density.

• dcompois2(): Conway-Maxwell-Poisson. Calculate density parameterized via the mean.

• pbinom(q = ad, size = ad, prob = ad, lower.tail = missing, log.p = missing): AD implementation ofpbinom

• pbinom(q = num, size = num, prob = num, lower.tail = missing, log.p = missing): Default method

• dmultinom(x = ad, size = ad., prob = ad, log = logical.): AD implementation ofdmultinom

• dmultinom(x = num, size = num., prob = num, log = logical.): Default method

• dmultinom(x = osa, size = ANY, prob = ANY, log = ANY): OSA implementation

• dmultinom(x = simref, size = ANY, prob = ANY, log = ANY): Simulation implementation. Modifiesx and returns zero.

• dmultinom(x = ANY, size = ANY, prob = ANY, log = ANY): Default implementation that checks for invalid usage.

• dcauchy(x = ad, location = ad., scale = ad., log = logical.): AD implementation ofdcauchy

• dcauchy(x = num, location = num., scale = num., log = logical.): Default method

• dcauchy(x = osa, location = ANY, scale = ANY, log = ANY): OSA implementation

• dcauchy(x = simref, location = ANY, scale = ANY, log = ANY): Simulation implementation. Modifiesx and returns zero.

• dgamma(x = ad, shape = ad, rate = ad, scale = missing, log = logical.): AD implementation ofdgamma

• pnbinom( q = ad, size = ad, prob = ad, mu = missing, lower.tail = missing, log.p = missing): AD implementation ofpnbinom

Examples

MakeTape( function(x) pnorm(x), x=numeric(5))$jacobian(1:5)

Interpolation

Description

Some interpolation methods are available to be used as part of 'RTMB' objective functions.

Usage

interpol1Dfun(z, xlim = c(1, length(z)), ...)interpol2Dfun(z, xlim = c(1, nrow(z)), ylim = c(1, ncol(z)), ...)## S4 method for signature 'ad,ad,ANY,missing'splinefun(x, y, method = c("fmm", "periodic", "natural"))## S4 method for signature 'ad,missing,ANY,missing'splinefun(x, method = c("fmm", "periodic", "natural"))

Arguments

Details

interpol1Dfun andinterpol2Dfun are kernel smoothers useful in the case where you need a 3rd ordersmooth representation of adata vector or matrix.A typical use case is when a high-resolution map needs to be accessed along a random effect trajectory.Both 1D and 2D cases accept an 'interpolation radius' parameter (default R=2) controlling the degree of smoothness. Note, that only the value R=1 will match the data exactly, while higher radius trades accuracy for smoothness. Note also that these smoothers do not attempt to extrapolate: The returned value will beNaN outside the valid range (xlim /ylim).

splinefun imitates the correspondingstats function. The AD implementation (in contrast tointerpol1Dfun) works for parameter dependent y-coordinates.

Value

function of x.

function of x and y.

Functions

• interpol1Dfun(): Construct a kernel smoothed representation of a vector.

• interpol2Dfun(): Construct a kernel smoothed representation of a matrix.

• splinefun(x = ad, y = ad, method = ANY, ties = missing): Construct a spline function.

• splinefun(x = ad, y = missing, method = ANY, ties = missing): Construct a spline function.

Examples

## ======= interpol1D## R=1 => exact match of observationsf <- interpol1Dfun(sin(1:10), R=1)layout(t(1:2))plot(sin(1:10))plot(f, 1, 10, add=TRUE)title("R=1")F <- MakeTape(f, 0)F3 <- F$jacfun()$jacfun()$jacfun()plot(Vectorize(F3), 1, 10)title("3rd derivative")## ======= interpol2D## R=1 => exact match of observationsf <- interpol2Dfun(volcano, xlim=c(0,1), ylim=c(0,1), R=1)f(0,0) == volcano[1,1]   ## Top-left cornerf(1,1) == volcano[87,61] ## Bottom-right corner## R=2 => trades accuracy for smoothnessf <- interpol2Dfun(volcano, xlim=c(0,1), ylim=c(0,1), R=2)f(0,0) - volcano[1,1]    ## Error Top-left cornerF <- MakeTape(function(x) f(x[1],x[2]), c(.5,.5))## ======= splinefunT <- MakeTape(function(x){   S <- splinefun(sin(x))   S(4:6)}, 1:10)

Multivariate Gaussian densities

Description

Multivariate Gaussian densities

Usage

dmvnorm(x, mu = 0, Sigma, log = FALSE, scale = 1)dgmrf(x, mu = 0, Q, log = FALSE, scale = 1)dautoreg(x, mu = 0, phi, log = FALSE, scale = 1)dseparable(...)unstructured(k)

Arguments

Details

Multivariate normal density evaluation is done usingdmvnorm(). This is meant for dense covariance matrices. Ifmany evaluations are needed for thesame covariance matrix please note that you can pass matrix arguments: Whenx is a matrix the density is applied to each row ofx and the return value will be a vector (length =nrow(x)) of densities.

The functiondgmrf() is essentially identical todmvnorm() with the only difference thatdgmrf() is specified via theprecision matrix (inverse covariance) assuming that this matrix issparse.

Autoregressive density evaluation is implemented for all orders viadautoreg() (including the simplest AR1).We note that this variant is for astationary,mean zero andvariance one process.FIXME: Provide parameterization via partial correlations.

Separable extension can be constructed for an unlimited number of inputs. Each input must be a function returning agaussianmean zerolog density. The output ofdseparable is anotherlog density which can be evaluated for array arguments. For exampledseparable(f1,f2,f3) takes as input a 3D arrayx.f1 acts in 1st array dimension ofx,f2 in 2nd dimension and so on. In addition tox, parametersmu andscale can be supplied - see below.

Value

Vector of densities.

Functions

• dmvnorm(): Multivariate normal distribution.OSA-residuals can be used for argumentx.

• dgmrf(): Multivariate normal distribution. OSA isnot implemented.

• dautoreg(): Gaussian stationary mean zero AR(k) density

• dseparable(): Separable extension of Gaussian log-densities

• unstructured(): Helper to generate an unstructured correlation matrix to use withdmvnorm

Scaling

All the densities accept ascale argument which replacesSCALE andVECSCALE functionality of TMB.Scaling is applied elementwise on the residualx-mu. This works as expected whenscale is ascalar or avector object of the same length asx.In addition,dmvnorm anddgmrf can be scaled by a vector of length equal to the covariance/precision dimension. In this case thescale parameter is recycled by row to meet the special row-wise vectorization of these densities.

Unstructured correlation

Replacement ofUNSTRUCTURED_CORR functionality of TMB. Constuct object usingus <- unstructured(k).Nowus has two methods:x <- us$parms() gives the parameter vector used as input to the objective function, andus$corr(x) turns the parameter vector into an unstructured correlation matrix.

Examples

func <- function(x, sd, parm, phi) {   ## IID N(0, sd^2)   f1 <- function(x)sum(dnorm(x, sd=sd, log=TRUE))   Sigma <- diag(2) + parm   ## MVNORM(0, Sigma)   f2 <- function(x)dmvnorm(x, Sigma=Sigma, log=TRUE)   ## AR(2) process   f3 <- function(x)dautoreg(x, phi=phi, log=TRUE)   ## Separable extension (implicit log=TRUE)   -dseparable(f1, f2, f3)(x)}parameters <- list(x = array(0, c(10, 2, 10)), sd=2, parm=1, phi=c(.9, -.2))obj <- MakeADFun(function(p)do.call(func, p), parameters, random="x")## Check that density integrates to 1obj$fn()## Check that integral is independent of the outer parametersobj$gr()## Check that we can simulate from this densitys <- obj$simulate()

Recursive quantile residuals

Description

OSA residuals are computed using the functiononeStepPredict. For this to work, you need to mark theobservation inside the objective function using theOBSfunction. Thereafter, residual calculation is as simple asoneStepPredict(obj). However, you probably want specify amethod to use.

Usage

oneStepPredict(  obj,  observation.name = names(obj$env$obs)[1],  data.term.indicator = "_RTMB_keep_",  ...)## S3 method for class 'osa'x[...]## S3 method for class 'osa'length(x)## S3 method for class 'osa'dim(x)## S3 method for class 'osa'is.array(x)## S3 method for class 'osa'is.matrix(x)

Arguments

Value

data.frame with standardized residuals; Same asoneStepPredict.

Functions

• oneStepPredict(): Calculate the residuals. See documentation ofTMB::oneStepPredict.

• [: Subset observations marked for OSA calculation.This function makes sure that when you subset an observation of class"osa" such asobs <- new("osa", x=advector(matrix(1:10,2)), keep = cbind(rep(TRUE,10),FALSE,FALSE))the 'keep' attribute will be adjusted accordinglyobs[,1:2]

• length(osa): Equivalent oflength

• dim(osa): Equivalent ofdim

• is.array(osa): Equivalent ofis.array

• is.matrix(osa): Equivalent ofis.matrix

Examples

set.seed(1)rw <- cumsum(.5*rnorm(20))obs <- rpois(20, lambda=exp(rw))func <- function(p) {  obs <- OBS(obs) ## Mark 'obs' for OSA calculation on request  ans <- 0  jump <- c(p$rw[1], diff(p$rw))  ans <- ans - sum(dnorm(jump, sd=p$sd, log=TRUE))  ans <- ans - sum(dpois(obs, lambda=exp(p$rw), log=TRUE))  ans}obj <- MakeADFun(func,                 parameters=list(rw=rep(0,20), sd=1),                 random="rw")nlminb(obj$par, obj$fn, obj$gr)res <- oneStepPredict(obj,                      method="oneStepGeneric",                      discrete=TRUE,                      range=c(0,Inf))$residual

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 theorder of likelihood accumulation in the user template - see details. Simulations can be obtained directly from the model object byobj$simulate() or used indirectly viacheckConsistency.In both cases a simulation is carried out from thefull model i.e. bothrandom effects andobservations are simulated. TheOBS function is used to mark which data items are observations - see examples.

Usage

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

Arguments

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 - seemethods(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 ofu[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 beNA).

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.

Value

An object with write access to store the simulation.

Functions

• simref(): Constructsimref

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

• length(simref): Equivalent oflength

• dim(simref): Equivalent ofdim

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

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

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

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

• [: Equivalent of[

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

• Ops(simref): Equivalent ofOps

• Math(simref): Equivalent ofMath

• t(simref): Equivalent oft

• diff(simref): Equivalent ofdiff

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

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 effectsfunc <- 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 workss <- simref(4)s2 <- 2 * s[1:2] + 1s2[] <- 7s ## 3 3 NA NA

Interface to TMB

Description

Interface to TMB

Usage

MakeADFun(  func,  parameters,  random = NULL,  profile = NULL,  integrate = NULL,  intern = FALSE,  map = list(),  ADreport = FALSE,  silent = FALSE,  ridge.correct = FALSE,  ...)sdreport(obj, ...)ADREPORT(x)REPORT(x)getAll(..., warn = TRUE)OBS(x)checkConsistency(obj, fast = TRUE, ...)

Arguments

Details

MakeADFun builds a TMB model object mostly compatible with theTMB package and with an almost identical interface.The main difference inRTMB is that the objective functionand the data is now given via a single argumentfunc. Becausefunc can be aclosure, there is no need for an explicit data argument toMakeADFun (see examples).

Value

TMB model object.

Functions

• MakeADFun(): Interface toMakeADFun.

• sdreport(): Interface tosdreport.

• ADREPORT(): Can be used inside the objective function to report quantities for which uncertainties will be calculated bysdreport.

• REPORT(): Can be used inside the objective function to report quantities via the model object usingobj$report().

• getAll(): Can be used to assign all parameter or data objects from a list inside the objective function.

• OBS(): Mark the observation to be used by eitheroneStepPredict or byobj$simulate.If your objective function is using an observationx, you simply needto runx <- OBS(x)inside the objective function.This will (1) allowoneStepPredict to change the class ofx to"osa" (OSA-residuals) or (2) allowobj$simulate to change the class ofx to"simref" (Simulation) on request.

• checkConsistency(): Interface tocheckConsistency.

Examples

## Objective with data from the user workspacedata(rivers)f <- function(p) { -sum(dnorm(rivers, p$mu, p$sd, log=TRUE)) }obj <- MakeADFun(f, list(mu=0, sd=1), silent=TRUE)opt <- nlminb(obj$par, obj$fn, obj$gr)sdreport(obj)## Same objective with an explicit data argumentf <- function(p, data) { -sum(dnorm(data, p$mu, p$sd, log=TRUE)) }cmb <- function(f, d) function(p) f(p, d) ## Helper to make closureobj <- MakeADFun(cmb(f, rivers), list(mu=0, sd=1), silent=TRUE)## 'REML trick'obj2 <- MakeADFun(cmb(f, rivers), list(mu=0, sd=1), random="mu", silent=TRUE)opt2 <- nlminb(obj2$par, obj2$fn, obj2$gr)sdreport(obj2) ## Compare with sd(rivers)## Single argument vector function with numeric 'parameters'fr <- function(x) {   ## Rosenbrock Banana function    x1 <- x[1]    x2 <- x[2]    100 * (x2 - x1 * x1)^2 + (1 - x1)^2}obj <- MakeADFun(fr, numeric(2), silent=TRUE)nlminb(c(-1.2, 1), obj$fn, obj$gr, obj$he)

The AD tape

Description

The AD tape as an R function

Usage

MakeTape(f, x)## S3 method for class 'Tape'x$name## S3 method for class 'Tape'print(x, ...)TapeConfig(  ...,  comparison = c("NA", "forbid", "tape", "allow"),  atomic = c("NA", "enable", "disable"),  vectorize = c("NA", "disable", "enable"))DataEval(f, x)GetTape(obj, name = c("ADFun", "ADGrad", "ADHess"), warn = TRUE)

Arguments

Details

A 'Tape' is a representation of a function that acceptsfixed size numeric input and returnsfixed size numeric output.The tape can be constructed usingF <- MakeTape(f, x) wheref is a standarddifferentiable R function (or more precisely: One using only functions that are documented to work for AD types).Having constructed a tape F, a number of methods are available:

Evaluation:

• Normal function evaluation 'F(x)' for numeric input.

• AD evaluation 'F(x)' as part of other tapes.

• Jacobian calculations using 'F$jacobian(x)'.

• Signal that the next forward pass should drop lazy evaluation and update the entire tapeF$force.update().

Transformation:

• Get new tape representing the Jacobian usingF$jacfun().

• Get new tape representing the sparse Jacobian usingF$jacfun(sparse=TRUE).

• Get new tape representing the Laplace approximation usingF$laplace(indices).

• Get new tape representing the Saddle Point approximation usingF$laplace(indices,SPA=TRUE).

• Get new tape representing the optimum (minimum) wrtindices byF$newton(indices).

• Get a 'shared pointer' representation of a tape usingF$atomic().

• Get tape of a single node byF$node(index) (mainly useful for derivative debugging).

• Reorder tape so selected inputs (indices) and its forward dependencies come as late as possible byF$reorder(indices).

Modification:

• Simplify internal representation of a tape usingF$simplify().

Extract tape information:

• Get internal parameter vector byF$par().

• Get computational graph byF$graph().

• Print the tape byF$print().

• Get internal arrays as adata.frame byF$data.frame().

• Find out where the tape spends its time byF$timer().

Value

Object of class"Tape".

Methods (by generic)

• $: Get a tape method.

• print(Tape): Print method

Functions

• MakeTape(): Generate a 'Tape' of an R function.

• TapeConfig(): Global configuration parameters of the tape (experts only!)comparison By default, AD comparison gives an error(comparison="forbid").This is the safe and recommended behaviour, because comparison is anon-differentiable operation. If you are building a tape thatrequires indicator functions e.g.f(x)*(x<0)+g(x)*(x>=0)then usecomparison="tape" to add the indicators to thetape. A final optioncomparison="allow" exists fortesting/illustration purposes. Do not use.

• DataEval(): Move a chunk of data from R to the tape by evaluating a normal R function (replaces TMB functionality 'DATA_UPDATE').

• GetTape(): Extract tapes from a model object created byMakeADFun.

Examples

F <- MakeTape(prod, numeric(3))show(F)F$print()H <- F$jacfun()$jacfun() ## Hessian tapeshow(H)#### Handy way to plot the graph of Fif (requireNamespace("igraph")) {   G <- igraph::graph_from_adjacency_matrix(F$graph())   plot(G, vertex.size=17, layout=igraph::layout_as_tree)}## Taped access of an element of 'rivers' datasetF <- MakeTape(function(i) DataEval( function(i) rivers[i] , i), 1 )F(1)F(2)## DATA_UPDATE example## - Pay attention to lazy evaluation!## - Use 'force.update()'mydat <- 1:4fetch <- function() { print("Getting 'mydat'"); .GlobalEnv$mydat }F <- MakeTape( function(x) x * sum(DataEval(fetch)), 1 )F$reorder() ## Re-order tape for faster executionF(1) ## No data updateF(2) ## No data updatemydat <- 5:8F(1) ## Still no data update !F$force.update() ## Signal that data has changedF(1) ## data updateF(2) ## No data update

Matrix exponential of sparse matrix multiplied by a vector.

Description

Calculatesexpm(A) %*% v using plain series summation. The number of terms is determined adaptively whenuniformization=TRUE.The uniformization method essentially pushes the spectrum of the operator inside a zero centered disc, within which a tight uniform error bound is available. This effectively reduces the number of terms needed to calculate the series to a given accuracy.IfA is a generator matrix (i.e.expm(A) is a probability matrix) and ifv is a probability vector, then the relative error of the result is bounded bytol.However, note that series summation may be unstable depending on the spectral radius of A. If you getNaN values, consider settingrescale_freq=1 for better stability (see details).

Usage

expAv(  A,  v,  transpose = FALSE,  uniformization = TRUE,  tol = 1e-08,  ...,  cache = A)

Arguments

Details

Additional supported arguments via... currently include:

• Nmax Integer (2e9 by default). Use no more than this number of terms even if the specified accuracy cannot be met. When usingexpAv as part of likelihood optimization, you can set a lower value to avoid long unnecessary computation when the optimizer tries extreme parameters. For example, if the spectral radius ofA cannot realistically exceed some known valuerhomax one can setNmax=qpois(tol, rhomax, lower.tail = FALSE).

• warn Logical (TRUE by default). Give warning if number of terms is truncated byNmax (autodiff code only).

• trace Logical (FALSE by default). Trace the number of terms when it adaptively changes (autodiff code only).

• rescale_freq Integer (50 by default) controlling how often to rescale for numerical stability. Set to a lower value for more frequent rescaling at the cost of longer computation time. The default value should suffice for a generator matrix of spectral radius up to at least1e6 (.Machine$double.xmax^(1/50)).

• rescale Logical; Set toFALSE to disable rescaling for higher speed. All other values are ignored.

Value

Vector (or matrix)

References

Grassmann, W. K. (1977). Transient solutions in Markovian queueing systems.Computers & Operations Research, 4(1), 47–53.

Sherlock, C. (2021). Direct statistical inference for finite Markov jump processes via the matrix exponential.Computational Statistics, 36(4), 2863–2887.

Examples

set.seed(1); A <- Matrix::rsparsematrix(5, 5, .5)expAv(A, 1:5) ## Matrix::expm(A) %*% 1:5F <- MakeTape(function(x) expAv(A*x, 1:5, trace=TRUE), 1)F(1)F(2) ## More terms needed => trigger retaping