另一个伪随机生成器干扰的伪随机序列

Question

我注意到，如果在生成伪随机序列时使用另一个伪随机数生成器，则会干扰种子序列。我的问题是是否有任何方法可以解决这个问题？你能否确保原始的种子序列继续下去？让我举个例子；以下是一个简单的for循环，它打印从正态分布中抽取的伪随机数：

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

以下是输出结果：

以下是输出结果：

[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

下面，如果迭代器等于五，我们介绍从均匀分布中产生的伪随机数。

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

以下是输出内容（在下面的输出中，星号表示从均匀分布中伪随机抽取的数）：

[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

我想要找到一个答案，即是否可以继续由set.seed(145)引入的原始种子序列，并因此获得以下输出：

[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

每一个输入都非常重要，特别是关于这个具体问题的文献参考。

编辑：

根据Rui Barradas的建议，我尝试自己实现一个函数，但没有成功。除了for循环内部的rnorm采样外，for循环内不应该有其他随机性，if语句中的采样应该由Rui的修复处理。但不幸的是，似乎有些东西正在干扰种子序列，因为下面的两个功能返回的不同，并且它们在如何绘制随机性（AR-1方程中通常的epsilon）的方式上有所不同，其他方面相同。

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)

我本来期望上述代码最终返回True，而不是告诉我平均相对差为1.186448的信息。

- Kristian Nielsen

Answer 1

系统变量.Random.seed存储了随机数生成器（RNG）的状态。从help(".Random.seed")中可以看到：

.Random.seed是一个整数向量，包含R中用于随机数生成的随机数生成器（RNG）状态。它可以保存和恢复，但不应由用户更改。

因此，以下内容有效。

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部分。

至于如何使此技巧在函数内工作，问题似乎是函数创建了它们自己的环境。而.Random.seed存在于.GlobalEnv中。因此，需要进行以下更改：使用.GlobalEnv$.Random.seed。

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

- Rui Barradas

---

谢谢，Rui！这正是我正在寻找的！我会立刻去阅读助手页面！ - Kristian Nielsen

你有时间给出你的看法吗，为什么这个技巧在函数内似乎不起作用？ - Kristian Nielsen

谢谢@RuiBarradas！不幸的是，在我的情况下似乎不起作用。我能否与您分享我的函数，因为我无法找到错误？ - Kristian Nielsen

@KristianNielsen 是的，您可以使用函数编辑问题，或者提出一个不同的问题，并附上指向此问题的链接。 - Rui Barradas

请查看我编辑后的问题@RuiBarradas - Kristian Nielsen

Answer 2

也许有更好的方法，但您可以预先计算随机值，然后在需要新值时引用该列表。以下将其转换为函数形式。您需要指定比最终所需的缓冲区大的缓冲区。这种方法的一个缺点是您需要提前指定随机函数和函数的参数。理论上，您可以使用反演变换抽样并从均匀分布中生成值来避免这个问题，但我将把它作为读者的练习...

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)
}

还有一个您可能会用到的示例...

> 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

- Dason

---

不可否认，这个函数的名称可以更好。不过在想到更好的名字之前，我需要喝更多的咖啡。 - Dason