tips for (new) RTMB users
- the objective function must be differentiable with respect to parameters (no
if(),abs(),round(),min(),max()depending on parameters) - have to implement exotic probability distributions yourself: anything that descends into C(++)/Fortran code will probably fail
- aliases to base functions made before call to
library(RTMB)might not work - use of
<-[(see here) etc.- specifically, if you use the
c()function, or if you use thediag<-function (which sets the diagonal of a matrix) or the[<-function (which assigns values within a matrix), you need to add e.g.ADoverload("[<-")to the beginning of your function
- specifically, if you use the
- You can only "grow" empty vectors if initialized with
numeric(0), notc()orNULL(i.e.x <- c(); x[1] <- 2throws an error, butx <- numeric(0); x[1] <- 2does work), see here. In any case it is better practice to preallocate instead of growing vectors, e.g.x <- numeric(1); x[1] <- 2(see chapter 2 of the R Inferno). - for matrix exponentials, you should use
Matrix::expm()rather thanexpm::expm() - RTMB is pickier than R about matrix types. You may need to use some combination of
drop()andas.matrix()to convert matrices with dimension 1 in some direction (orMatrixmatrices) back to vectors [[-indexing may be much faster than[-indexing: see here (and later messages in that thread)- if you use
cat()orprint()to print out numeric values, the results may not make sense (you'll see a printout of RTMB's internal representation of autodiff-augmented numbers ...)
- RTMB uses
%*%(as in base R), not*(as in C++) for matrix/matrix and matrix/vector multiplication
- as in C++ code, you probably don't have to worry as much about non-vectorized code (e.g.
forloops) being very slow relative to vectorized codes: seevectorization_comp.Rfor the code behind this example:
Unit: microseconds
expr min lq mean median uq max neval cld
vect 2.8 3.2 4.3 3.5 3.9 282 10000 a
vect_RTMB 12.7 13.4 15.0 13.8 14.8 135 10000 b
forloop 223.0 240.1 245.0 244.2 248.7 382 10000 c
forloop_RTMB 12.6 13.5 15.7 13.9 14.8 3835 10000 b
- data handling (see here, here) (and very similar arguments from 2004 about MLE fitting machinery taking a
dataargument) - have to handle prediction, tests, diagnostics, etc. etc. yourself (although the
broom.mixedpackage does have a rudimentarytidymethod for TMB objects (e.g.class(TMBobj) <- c("TMB", class(TMBobj)); broom.mixed::tidy(TMBobj)) - if you do something clever where you define your objective function in a different environment from where you call
MakeADFun, you can useassign(..., environment(objective_function))to make sure that the objective function can see any objects it needs to know about ...