Nothing
## ----setup, include=FALSE, message=FALSE--------------------------------------
##library(knitr)
knitr::opts_chunk$set(fig.align='center')
library(RTMB)
set.seed(1)
formals(MakeADFun)$silent <- TRUE
## Utilites for graph visulization
library(igraph)
## To tweak the vertex names
iformat <- function(n, to=c("superscript", "subscript")) {
to <- match.arg(to)
digits <- strsplit(as.character(n),"")[[1]]
if (to=="superscript")
x <- paste0("0x", c("2070", "00B9", "00B2", "00B3", paste0("207",4:9)))
else
x <- paste0("208", 0:9)
names(x) <- 0:9
intToUtf8(x[digits])
}
addindex <- function(x, to) {
paste0(x, sapply(seq_along(x), iformat, to))
}
showGraph <- function(F) {
g <- F$graph()
colnames(g)[colnames(g)=="Inv"] <- "X"
colnames(g)[colnames(g)=="Dep"] <- "Y"
colnames(g) <- addindex(colnames(g), "sup")
G <- graph_from_adjacency_matrix(g)
oldpar <- par(mar=c(0,0,0,0),oma=c(0,0,0,0))
on.exit(par(oldpar))
##plot(G, vertex.size=17, layout=layout_as_tree)
plot(G, vertex.size=23, layout=layout_as_tree)
}
## -----------------------------------------------------------------------------
f <- function(x) exp( x[1] + 1.23 * x[2] )
## -----------------------------------------------------------------------------
F <- MakeTape(f, numeric(2))
F
## -----------------------------------------------------------------------------
F$methods()
## -----------------------------------------------------------------------------
F$print()
## -----------------------------------------------------------------------------
F(3:4)
## -----------------------------------------------------------------------------
F$print()
## -----------------------------------------------------------------------------
F$jacobian(3:4)
## -----------------------------------------------------------------------------
F$print()
## -----------------------------------------------------------------------------
F$graph()
## ----fig.cap="Operator graph of test function"--------------------------------
showGraph(F)
## -----------------------------------------------------------------------------
G <- MakeTape(function(x) c( F(x) , F(x) ), numeric(2))
## ----fig.cap="Operator graph of test function evaluated twice"----------------
showGraph(G)
## -----------------------------------------------------------------------------
F <- F$atomic()
## ----fig.cap="Operator graph of *atomic* test function evaluated twice"-------
G <- MakeTape(function(x) c( F(x) , F(x) ), numeric(2))
showGraph(G)
## -----------------------------------------------------------------------------
F <- MakeTape(function(x) {
y1 <- sin(x[1] + x[2])
y2 <- sin(x[1] + x[2])
y3 <- cos(x[1] + x[2])
y1+y2
}, numeric(2))
## ----fig.cap="Tape of function"-----------------------------------------------
showGraph(F)
## -----------------------------------------------------------------------------
F$simplify("eliminate")
## ----fig.cap="Tape of function after eliminate"-------------------------------
showGraph(F)
## -----------------------------------------------------------------------------
F$simplify("optimize")
## ----fig.cap="Tape of function after optimize"--------------------------------
showGraph(F)
## -----------------------------------------------------------------------------
D <- MakeTape(diff, numeric(5))
D
## -----------------------------------------------------------------------------
J <- D$jacfun()
J
## -----------------------------------------------------------------------------
J(1:5)
## -----------------------------------------------------------------------------
Js <- D$jacfun(sparse=TRUE)
Js
## -----------------------------------------------------------------------------
Js(1:5)
## -----------------------------------------------------------------------------
f <- function(X) X %*% X
## ----fig.cap="Plain matrix multiply: Expands all operations"------------------
TapeConfig(atomic="disable")
F <- MakeTape(f, matrix(0, 2, 2))
showGraph(F)
## ----fig.cap="Atomic matrix multiply: Collapses to a single operation. The constants are the matrix dimensions which are represented as additional inputs."----
TapeConfig(atomic="enable")
F <- MakeTape(f, matrix(0, 2, 2))
showGraph(F)
## -----------------------------------------------------------------------------
F <- MakeTape(function(x) {
a <- sin(x)
G <- MakeTape(function(y) {
a * y
}, numeric(1))
DG <- G$jacfun()
DG(x * x)
}, numeric(1))
## -----------------------------------------------------------------------------
## Negative log of the integrand
f <- MakeTape(function(x)x[1]^2+x[2]^2, numeric(2))
show(f)
## -----------------------------------------------------------------------------
## Integrate x[2] out and return the negative log of the result
F <- f$laplace(2)
show(F)
## -----------------------------------------------------------------------------
F(3)
-log(integrate(function(x2)exp(-3^2-x2^2), -Inf, Inf)$value)
## -----------------------------------------------------------------------------
## minimize wrt. x[2] and return optimum as function of x[1]
S <- f$newton(2)
show(S)
## -----------------------------------------------------------------------------
S(3)
## -----------------------------------------------------------------------------
rdblpois <- function(n, muX=5, muN=10) {
replicate(n, sum(rpois(rpois(1, muN), muX)))
}
set.seed(1)
x <- rdblpois(100)
## -----------------------------------------------------------------------------
Kpois <- function(t, mu) mu * (exp(t) - 1)
## -----------------------------------------------------------------------------
Kpois2 <- function(s, muX, muN) Kpois(Kpois(s, muX), muN)
## -----------------------------------------------------------------------------
nldens <- function(obs, muX, muN) {
## 1
F <- MakeTape(function(s) {
K <- Kpois2(s, muX, muN) ## CGF
K <- K - s * obs ## SPA adjustment
sum(K)
}, rep(1, length(obs)))
## 2
L <- F$laplace(1:length(obs), SPA=TRUE)
## 3
L(numeric(0))
}
## -----------------------------------------------------------------------------
obj <- MakeADFun(function(p) nldens(x, p$muX, p$muN), list(muX=1, muN=1), silent=TRUE)
opt <- nlminb(obj$par, obj$fn, obj$gr)
sdreport(obj)
## -----------------------------------------------------------------------------
f <- function(x) {
## Get real/imag part
xreal <- x[1]
ximag <- x[2]
## Construct AD complex number
z <- as.complex(xreal) + 1i * as.complex(ximag)
## Return real numbers
Mod(exp(z))
}
## -----------------------------------------------------------------------------
F <- MakeTape(f, numeric(2))
F
## -----------------------------------------------------------------------------
## Using R complex arithmetic
f(1:2)
## -----------------------------------------------------------------------------
## Using the tape representation
F(1:2)
## -----------------------------------------------------------------------------
F$print()
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.