Pseudo-random sequence interfered by another pseudo-random generator

Question

I have noticed that if one uses another pseudo-random number generator when generating a pseudo-random sequence, the seed sequence is interfered. My question is if there is anything to do about it? Can you somehow ensure original seed sequence is continued? Let me show an example;

A simple for-loop which prints a pseudo-random number drawn from the normal distribution:

set.seed(145)
for (i in 1:10){
  print(rnorm(1,0,1))
}

Which gives the following output:

[1] 0.6869129
[1] 1.066363
[1] 0.5367006
[1] 1.906029
[1] 1.06316
[1] 1.370344
[1] 0.5277918
[1] 0.4030967
[1] 1.167752
[1] 0.7926794

Next, we introduce a pseudo-random draw from the uniform distribution if the iterator is equal to five.

set.seed(145)
for (i in 1:10){
  print(rnorm(1,0,1))
  if (i == 5){
    print(runif(1,0,1))
  }
}

Which gives the following output (in the following output, the star marks the pseudo-random draw from the uniform distribution):

[1] 0.6869129
[1] 1.066363
[1] 0.5367006
[1] 1.906029
[1] 1.06316
[1] 0.9147102*
[1] -1.508828
[1] -0.03101992 
[1] -1.091504
[1] 0.2442405
[1] -0.6103299

What I try to seek an answer on, is whether or not it is possible to continue the original seed sequence introduced by set.seed(145), and thereby get the following output:

[1] 0.6869129
[1] 1.066363
[1] 0.5367006
[1] 1.906029
[1] 1.06316
[1] 0.9147102*
[1] 1.370344
[1] 0.5277918
[1] 0.4030967
[1] 1.167752
[1] 0.7926794

Every input is highly appreciated, especially some references to literature on this specific issue.

EDIT:

Based on the input of Rui Barradas I tried to implement it in a function of my own but without luck. Besides the rnorm sampling within each iteration of the for-loop, there shouldn't be any other randomness in the for-loop expect for the sampling in the if-statement which should be handled by the fix of Rui. But unfortunately there seems to be something interfering the seed sequence, as the two functions below does not return the same, and they are equal with the exception of how the randomness (epsilon normally in AR-1 equation) is drawn.

tt <- rnorm(500,0,1)*10 

test1 <- function(y, x0=1, n,qsigma = 3, alpha = 5, beta = 20, limit = 0.30){
  t <- length(y)
  gama <- (alpha + beta)/2
  x <- matrix(0,n,t)
  x[, 1] <- rep(x0,n)
  for(s in 2:t) {
    x[, s] <-pmax(alpha*(x[,s-1]<=gama) +beta*(x[,s-1]>gama)+rnorm(n,0,qsigma),1)
    if (s==250) {
      current <- .GlobalEnv$.Random.seed
      resamp <- sample(n, n, replace = TRUE)
      x[,s] <- x[resamp,s]
      .GlobalEnv$.Random.seed <- current
      }
  }
  list(x = x)
}

test3 <- function(y, x0=1, n,qsigma = 3, alpha = 5, beta = 20, limit = 0.30) {
  t <- length(y)
  gama <- (alpha + beta)/2
  x <- matrix(0,n,t)
  x[, 1] <- rep(x0,n)
  e_4 <- matrix(rnorm(n * (t), 0, qsigma),n, (t))

  for(s in 2:t) {
    x[, s] <-pmax(alpha*(x[,s-1]<=gama) +beta*(x[,s-1]>gama)+e_4[,(s-1)],1)
    if (s==250) {resamp <-sample(n, n, replace = TRUE)
      x[,s] <- x[resamp,s]
    }
  }
  list(x = x, pp = e_4)
}

set.seed(123)
dej11 <- test3(y = tt, n = 5000)$x
set.seed(123)
dej21 <- test1(y = tt, n = 5000)$x
all.equal(dej11,dej21)

I did expect the above to in the end return True instead of a message telling me that the mean relative difference is 1.186448.

Answer 1

The system variable .Random.seed stores the state of the rng. From help(".Random.seed"):

.Random.seed is an integer vector, containing the random number generator (RNG) state for random number generation in R. It can be saved and restored, but should not be altered by the user.

So the following works.

set.seed(145)
for (i in 1:10){
  print(rnorm(1,0,1))
  if (i == 5){
    current <- .Random.seed
    print(runif(1,0,1))
    .Random.seed <- current
  }
}

Note that you should read that help page carefully, in particular section Note.

As for how to make this trick work inside a function, the problem seems to be that functions create their own environments. And the .Random.seed exists in the .GlobalEnv. So the following change is needed: use .GlobalEnv$.Random.seed instead.

set.seed(145)

f <- function() {
    for (i in 1:10) {
        print(rnorm(1, 0, 1))
        if (i == 5) {
            current <- .GlobalEnv$.Random.seed
            print(runif(1, 0, 1))
            .GlobalEnv$.Random.seed <- current
        }
    }
}

f()
#[1] 0.6869129
#[1] 1.066363
#[1] 0.5367006
#[1] 1.906029
#[1] 1.06316
#[1] 0.9147102
#[1] 1.370344
#[1] 0.5277918
#[1] 0.4030967
#[1] 1.167752
#[1] 0.7926794

Answer 2

There are probably better ways but you could just pre-compute your random values and then refer to that list when you need new values. The following will put that into a function form. You will want to specify a buffer larger than what you ultimately will need. One disadvantage of this approach is that you do need to specify the random function and the parameters for the function ahead of time. Theoretically you could use inverse transform sampling and just generate values from a uniform distribution to get around this but I'll leave that as an exercise for the reader...

random_seed_fixed <- function(rfun, seed, buffer = 1000000, ...){
  set.seed(seed)
  values <- rfun(buffer, ...)
  next_index <- 1

  out <- function(n){
    new_index <- next_index + n
    # Give an error if we're going to exceed the bounds of our values
    stopifnot(new_index < buffer)

    id <- seq(next_index, new_index - 1, by = 1)
    next_index <<- new_index
    ans <- values[id]
    return(ans)
  }

  return(out)
}

And an example of how you might use it...

> my_rnorm <- random_seed_fixed(rnorm, seed = 642, mean = 17, sd = 2.3)
> 
> my_rnorm(5)
[1] 18.53370 16.16721 15.43144 16.67967 18.27675
> my_rnorm(5)
[1] 19.26933 17.50994 18.90019 14.80153 18.18837
> 
> my_rnorm <- random_seed_fixed(rnorm, seed = 642, mean = 17, sd = 2.3)
> my_rnorm(5) # matches the previous first call of my_rnorm(5)
[1] 18.53370 16.16721 15.43144 16.67967 18.27675
> rnorm(1, 0, 1)
[1] 2.515765
> my_rnorm(5) # Still matches the previous second call of my_rnorm(5)
[1] 19.26933 17.50994 18.90019 14.80153 18.18837