Random numbers

Overview

ast2ast has no global seed. Random numbers come from engine and distribution objects that live inside the translated function. The numbers are the same as those of R after set.seed(seed) with the default RNGkind(). normal and dnorm can differ from R in the last bit when the compiler uses fused multiply-add (FMA) and R was built without it.

Drawing numbers

f <- function(seed) {
  argtypes(seed |> type(int))
  e <- mersenne_twister(seed)
  d <- normal(e, mean = 2, sd = 0.5)
  draw(d, 5L)
}
fcpp <- ast2ast::translate(f)
fcpp(1L)

set.seed(1)
rnorm(5, 2, 0.5) # same numbers

Every draw() advances the engine of the distribution, so two calls give different numbers, exactly as two calls of rnorm.

mean, sd, min, max and seed can be read with $ (d$mean, e$seed) but not written.

Sampling

sample and sample_int take the engine as first argument. All three variants of R are supported: without replacement, with replacement and with prob.

f <- function(x, p) {
  argtypes(x |> type(vec(double)), p |> type(vec(double)))
  e <- mersenne_twister(6L)
  sample(e, x, 10L, TRUE, p)
}
x <- c(1.5, 2.5, 3.5, 4.5, 5.5)
p <- c(1, 1, 2, 2, 1)
ast2ast::translate(f)(x, p)

set.seed(6)
sample(x, 10, TRUE, p) # identical

As in R, a single number x >= 1 means sample(1:x).

Densities

dnorm and dunif are recycled as in R; the result is a scalar if x and the parameters are scalars. They work in every translation; dnorm has derivatives w.r.t. x, mean and sd, dunif w.r.t. min and max. This makes them usable in log-likelihoods for lbfgsb, hessian or HMC.

nll <- function(theta, y) {
  argtypes(theta |> type(vec(double)), y |> type(vec(double)))
  -sum(dnorm(y, theta[[1L]], exp(theta[[2L]]), log = TRUE))
}
fcpp <- ast2ast::translate(nll, derivative = "reverse")

set.seed(1)
y <- rnorm(50, mean = 3, sd = 2)
theta <- c(3, log(2))
fcpp(theta, y)
-sum(dnorm(y, 3, 2, log = TRUE)) # same value

Each distribution owns its engine

normal(e) and uniform(e) store a copy of e. Consequences:

• Two distributions created from the same engine produce the same stream.

  f <- function() {
    e <- mersenne_twister(3L)
    a <- normal(e)
    b <- normal(e)
    c(draw(a, 3L), draw(b, 3L)) # rep(rnorm(3), 2) after set.seed(3)
  }

  Use different seeds for independent streams.

• Drawing from a does not advance e, and sample(e, ...) does not advance a.

• Mixing distributions is therefore not identical to R. In R, rnorm(2); runif(2) share one stream; in ast2ast each distribution has its own.

Passing to inner functions

Pass engines and distributions with ref(). Passed by value, the inner function draws from a copy, and the caller later draws the same numbers again. ast2ast warns at translation time in this case.

f <- function() {
  g <- function(d) {
    argtypes(d |> type(normal) |> ref())
    returntype(vec(double))
    draw(d, 2L)
  }
  d <- normal(mersenne_twister(9L))
  c(g(d), g(d)) # rnorm(4) after set.seed(9)
}

A const() parameter is an error, since drawing changes the object.

Order of evaluation

C++ does not specify the order in which function arguments are evaluated (GCC evaluates right to left); R evaluates left to right. With several draws in one expression the result can therefore differ from R:

f <- function() {
  d <- normal(mersenne_twister(1L))
  c(draw(d, 2L), draw(d, 2L)) # the two halves may be swapped
}

Assign each draw to a variable first:

f <- function() {
  d <- normal(mersenne_twister(1L))
  a <- draw(d, 2L)
  b <- draw(d, 2L)
  c(a, b) # rnorm(4) after set.seed(1)
}

Limitations

• Engines and distributions cannot be passed from R or returned to R.
• The type names mersenne_twister, normal and uniform are reserved; new_type() cannot reuse them.
• Only the Mersenne-Twister engine and the normal and uniform distribution are available. No p*/q* functions yet.
• sample without replacement and without prob differs from R for x > 1e7 and size <= x / 2 (R then uses hashing).