Reproducible Randomness Is More Than Calling set.seed()

Calling set.seed() is usually the first thing we learn about reproducible simulation in R.

It is necessary often enough that it becomes almost synonymous with reproducibility:

1
2
set.seed(42)
rnorm(5)

Run the code again with the same R version, RNG configuration and seed, and the same pseudo-random sequence is generated.

That is useful.

It is not yet a good random-number contract for a statistical function.

Consider a package function that runs a bootstrap internally:

1
2
3
4
5
6
7
8
9
10
bootstrap_mean <- function(x, B = 999L, seed = 1L) {
  set.seed(seed)

  estimates <- replicate(
    B,
    mean(sample(x, replace = TRUE))
  )

  quantile(estimates, c(0.025, 0.975))
}

The function is reproducible.

It also silently replaces the caller's random-number stream.

That second fact is easy to miss because the function returns the right object. The damage appears later, in code that has nothing visibly to do with the function call.

For statistical software, I want a stronger guarantee:

The function should control its own randomness when explicitly asked to do so, without controlling the caller's future randomness as a side effect.

The Hidden Global Variable

R's ordinary pseudo-random-number generator is stateful.

After random numbers have been generated, its state is represented by the object

1
.Random.seed

in the global environment.

A call such as

1
runif(1)

uses the current state and replaces it with the next state in the stream.

Schematically,

where (S_t) is the generator state and (X_t) the generated value.

set.seed(42) does something different. It replaces the current state with the state associated with seed 42:

That is exactly what we want at the top of a script when we intentionally define the random experiment.

Inside a reusable function it is a global side effect.

A Reproducible Function Can Break a Reproducible Script

Here is a smaller example.

1
2
3
4
naive_random_mean <- function(seed) {
  set.seed(seed)
  mean(rnorm(100))
}

Now compare a caller's intended random stream with the stream after that function is inserted between two draws.

1
2
3
4
5
6
7
8
9
10
11
12
set.seed(100)
expected <- c(runif(1), runif(1))

set.seed(100)
first <- runif(1)

invisible(naive_random_mean(42))

second <- runif(1)

identical(c(first, second), expected)
# FALSE

The caller seeded the experiment correctly.

The function changed it anyway.

The failure is compositional. Each piece of code looks reproducible in isolation, but combining them changes the experiment.

This is why I do not think the contract

same seed, same result

is sufficient for a statistical package.

A reusable function also needs to say what happens to the surrounding random stream.

The Four Invariants I Want

For a statistical function with an optional seed argument, I would test four behaviours.

1. Explicit seed, existing caller state

If the caller already has an RNG state and supplies a function-local seed, the caller state should be byte-for-byte identical after the call:

2. Explicit seed, no existing caller state

If .Random.seed did not exist before the call, it should not suddenly exist afterwards.

The function should not leave random-generator state behind merely because its internal implementation used simulation.

3. No explicit seed

If seed = NULL, the function should normally consume the caller's current random stream.

Then

when the function actually draws randomness.

This preserves the usual R semantics. An unseeded stochastic call belongs to the caller's experiment.

4. Same explicit seed, same result

Two calls with the same inputs and explicit seed should reproduce the same stochastic result.

This is the familiar part, but it belongs together with the other three.

The full contract is therefore:

A Local Seed Wrapper in Base R

A small wrapper can implement that contract.

I prefer a callback here because it makes the evaluation boundary explicit.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
with_local_seed <- function(seed, fn) {
  stopifnot(is.function(fn))

  if (is.null(seed)) {
    return(fn())
  }

  if (
    length(seed) != 1L ||
    is.na(seed) ||
    seed != as.integer(seed)
  ) {
    stop("`seed` must be NULL or a single integer.", call. = FALSE)
  }

  old_kind <- RNGkind()
  had_seed <- exists(
    ".Random.seed",
    envir = .GlobalEnv,
    inherits = FALSE
  )

  if (had_seed) {
    old_seed <- get(
      ".Random.seed",
      envir = .GlobalEnv,
      inherits = FALSE
    )
  }

  on.exit({
    # Restore the generator configuration before restoring its exact state.
    do.call(RNGkind, as.list(old_kind))

    if (had_seed) {
      assign(
        ".Random.seed",
        old_seed,
        envir = .GlobalEnv
      )
    } else if (
      exists(
        ".Random.seed",
        envir = .GlobalEnv,
        inherits = FALSE
      )
    ) {
      rm(".Random.seed", envir = .GlobalEnv)
    }
  }, add = TRUE)

  set.seed(as.integer(seed))
  fn()
}

The on.exit() is the important part.

The state must be restored not only after a successful result but also after an error.

The wrapper also records RNGkind(). A function that only calls set.seed() does not change the RNG kind, so this is slightly more defensive than strictly necessary for the minimal example. It becomes necessary if code inside the boundary changes the generator, normal generator or sampling algorithm.

A Better Monte Carlo Function

Now the stochastic function can make its seed semantics explicit.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
random_mean <- function(n = 100L, seed = NULL) {
  if (
    length(n) != 1L ||
    is.na(n) ||
    n < 1L ||
    n != as.integer(n)
  ) {
    stop("`n` must be a positive integer.", call. = FALSE)
  }

  with_local_seed(seed, function() {
    mean(rnorm(as.integer(n)))
  })
}

A seeded call is now local:

1
2
3
4
5
6
7
8
9
10
11
12
set.seed(100)
expected <- c(runif(1), runif(1))

set.seed(100)
first <- runif(1)

invisible(random_mean(seed = 42L))

second <- runif(1)

identical(c(first, second), expected)
# TRUE

The function gets deterministic internal randomness while the caller keeps the stream it already had.

That is much closer to how I expect a library function to behave.

Reproducibility Without State Leakage

The simplest regression checks are direct.

Existing state is restored exactly

1
2
3
4
5
6
7
8
9
set.seed(123L)
before <- .Random.seed

invisible(random_mean(seed = 7L))

after <- .Random.seed

identical(after, before)
# TRUE

I would use identical(), not a numeric tolerance.

The RNG state is discrete program state. If the contract says it is restored, approximate equality is not the claim.

No state is left behind

This branch deserves its own test.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
if (
  exists(
    ".Random.seed",
    envir = .GlobalEnv,
    inherits = FALSE
  )
) {
  rm(".Random.seed", envir = .GlobalEnv)
}

invisible(random_mean(seed = 7L))

exists(
  ".Random.seed",
  envir = .GlobalEnv,
  inherits = FALSE
)
# FALSE

This case is easy to forget because most interactive R sessions already have a .Random.seed by the time tests are inspected manually.

But it is a different branch of the state-management logic.

If it is wrong, a function that claims to be locally seeded can leak a new global state into a previously clean session.

No seed means normal stream consumption

The opposite behaviour matters too.

1
2
3
4
5
6
7
8
9
set.seed(321L)
before <- .Random.seed

invisible(random_mean(seed = NULL))

after <- .Random.seed

identical(after, before)
# FALSE

A wrapper that always restores the caller state would be wrong here.

It would make an apparently random unseeded function repeat the same draws whenever called from the same point in the stream.

Locality should be opt-in through the explicit seed, not imposed on ordinary stochastic execution.

Same explicit seed means exact replay

1
2
3
4
5
first <- random_mean(n = 1_000L, seed = 42L)
second <- random_mean(n = 1_000L, seed = 42L)

identical(first, second)
# TRUE

This is the usual reproducibility assertion.

It is now one quarter of the contract rather than the whole contract.

Error Paths Are Part of the RNG Contract

State restoration that works only when the simulation succeeds is not state restoration.

Suppose the random work throws an error after drawing some values:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
set.seed(12L)
before <- .Random.seed

try(
  with_local_seed(99L, function() {
    rnorm(10)
    stop("simulation failed")
  }),
  silent = TRUE
)

after <- .Random.seed

identical(after, before)
# TRUE

This is why on.exit() is a better primitive than manually restoring the state on the final line of a function.

Statistical code has many ways to exit early:

• failed optimisers,
• singular matrix decompositions,
• non-finite simulated statistics,
• interrupted loops,
• explicit validation errors after stochastic preprocessing.

The global RNG state should not depend on which exit path happened to fire.

Why set.seed() at the Top of Every Function Is a Bad Pattern

It is sometimes defended as a way to make package behaviour stable:

1
2
3
4
my_bootstrap <- function(x) {
  set.seed(123)
  # ...
}

This creates two problems.

First, every call gets the same random sequence whether or not the caller asked for it.

Second, every call resets the global stream.

A stochastic algorithm with hidden fixed seeding is not really random from the caller's perspective. Repeated calls can reuse the same resamples, Monte Carlo draws or initialisations in ways that are not visible in the API.

A seed is part of the experimental design. It should not be an undocumented constant buried inside a function.

A better interface is

1
my_bootstrap(x, seed = NULL)

with documented behaviour for both NULL and an explicit integer.

Why Saving Only the Integer Seed Is Not Enough

Another tempting pattern is:

1
2
3
set.seed(seed)
# do work
set.seed(old_seed)

There is no old_seed integer that generally represents the caller's current point in the random stream.

The full state is .Random.seed, an integer vector encoding both generator information and its current state.

The caller may have executed thousands of random draws since the last explicit set.seed() call.

What must be restored is the state itself:

1
old_seed <- .Random.seed

not merely some earlier seed value.

This is an important distinction:

The seed initializes a stream. The state identifies where we currently are in it.

The Test Can Be Wrong Too

R adds one more subtlety: arguments are evaluated lazily.

That matters when the object used by an RNG test is itself created by code that changes the seed.

Consider this helper:

1
2
3
4
make_input <- function() {
  set.seed(1L)
  rnorm(20)
}

and a function that accepts an input as a lazy argument:

1
2
3
make_runner <- function(x) {
  function() mean(x)
}

This looks harmless:

1
runner <- make_runner(make_input())

But make_input() may still be an unevaluated promise. It can be forced later, when runner() first needs x.

An RNG-state test can therefore do this accidentally:

1. construct a lazy fixture that will call set.seed() later,
2. record .Random.seed,
3. call the function under test,
4. force the fixture inside that call,
5. observe that the RNG state changed,
6. blame the function under test.

The test has measured its own fixture.

The fix is simple when a helper is meant to capture a fully constructed object:

1
2
3
4
5
make_runner <- function(x) {
  force(x)

  function() mean(x)
}

The broader lesson is better:

Testing random-state side effects requires thinking about evaluation order, not just values.

Test the State, Not Only the Result

Many stochastic tests look like this:

1
2
3
4
5
6
7
set.seed(42)
a <- my_function()

set.seed(42)
b <- my_function()

expect_equal(a, b)

That verifies replay.

It says nothing about whether my_function() damaged the surrounding stream.

For package-quality testing, I would add explicit state assertions.

Using testthat, the four core cases look like this:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
test_that("explicit seed restores existing RNG state", {
  set.seed(123L)
  before <- .Random.seed

  invisible(random_mean(seed = 7L))

  expect_identical(.Random.seed, before)
})


test_that("explicit seed leaves no state when none existed", {
  if (exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE)) {
    rm(".Random.seed", envir = .GlobalEnv)
  }

  invisible(random_mean(seed = 7L))

  expect_false(
    exists(".Random.seed", envir = .GlobalEnv, inherits = FALSE)
  )
})


test_that("NULL seed advances caller stream", {
  set.seed(321L)
  before <- .Random.seed

  invisible(random_mean(seed = NULL))

  expect_false(identical(.Random.seed, before))
})


test_that("same explicit seed reproduces result", {
  expect_identical(
    random_mean(1_000L, seed = 42L),
    random_mean(1_000L, seed = 42L)
  )
})

I would add an error-path test as a fifth assertion whenever the local-seed helper is part of package infrastructure.

Monte Carlo Replicates Are Part of Reproducibility Too

A seed alone does not define a Monte Carlo procedure.

Suppose a p-value is estimated as

where (K) is the number of simulated statistics at least as extreme as the observed statistic.

Then both

and

are part of the computational experiment.

Changing (B) changes the set of possible p-values:

For a bootstrap confidence interval, the replicate count similarly affects Monte Carlo variability in the estimated quantiles.

So reproducibility metadata should retain at least:

• the seed,
• the number of replicates,
• algorithm settings that affect random draws,
• the software version when exact replay matters.

set.seed(42) by itself is not a complete computational specification.

Seeds Should Be Inputs, Not Provenance Lost in a Script

If stochastic output matters scientifically, I prefer the seed to appear explicitly in the result or experiment record.

For example:

1
2
3
4
5
6
result <- list(
  estimate = estimate,
  interval = interval,
  B = B,
  seed = seed
)

That does not mean every end-user print method needs to display the seed prominently.

It means the experiment can be reconstructed without guessing which line in a script happened to initialize the generator three pages earlier.

A random seed is small provenance with unusually high value.

Parallel Randomness Is a Separate Problem

Everything above concerns one ordinary R process.

Parallel simulation adds another layer.

Naively sending the same seed to several workers can create identical streams. Letting workers inherit state can make results depend on scheduling or process-launch details.

For parallel Monte Carlo, independent reproducible streams should be designed explicitly. In base R that usually points toward the L'Ecuyer-CMRG generator family and its stream/substream machinery. Higher-level parallel frameworks often provide their own reproducible seeding contracts.

The principle remains the same:

A local-seed wrapper solves scope inside one process. It is not a substitute for a parallel RNG design.

Randomness Is Part of the API

When a statistical function contains a bootstrap, permutation test, Monte Carlo integral, random initialization or simulation step, its RNG behaviour is part of its public interface whether the documentation acknowledges it or not.

Users need to know:

• whether repeated calls are stochastic,
• whether an explicit seed can reproduce a result,
• whether that seed is local to the function,
• whether an unseeded call consumes the current caller stream,
• whether parallel execution changes the guarantee.

Those are API semantics, not implementation trivia.

A function that returns numerically correct output but invisibly resets the global RNG can make the next analysis numerically wrong relative to the experiment the user intended to run.

The Practical Rule

For scripts, notebooks and one-off experiments, this remains perfectly sensible:

1
set.seed(42)

Set the seed at the experiment boundary and make the whole sequence reproducible.

For reusable statistical functions, I would use a different rule:

That means snapshotting the caller's state, running the stochastic computation under the requested seed, and restoring the state even if the computation fails.

If no seed is supplied, let the function participate normally in the caller's random stream.

That gives us the property we actually want:

That is a stronger form of reproducibility than merely calling set.seed().