R中rbinom（prob = 0.5）的不稳定种子行为

Posted 2023-02-25

【中文标题】R中rbinom（prob = 0.5）的不稳定种子行为

【英文标题】：Erratic seed behavior with rbinom(prob=0.5) in R

【发布时间】：2013-09-25 07:12:48

【问题描述】：

当R 与rbinom() 一起使用prob=0.5 时，我发现了我认为不稳定的行为（但我希望有一个简单的解释）。总体思路：对我来说，如果我设置种子，运行一次rbinom()（即执行一个随机过程），尽管prob 设置为什么值，随机 种子应该改变一个增量。然后，如果我再次将种子设置为相同的值，并运行另一个随机过程（例如再次rbinom()，但可能使用不同的值prob），种子应该再次更改为相同的值对于前面的单个随机过程。

我发现R 正是这样做的，只要我将rbinom() 与任何prob!=0.5 一起使用。这是一个例子：

比较种子向量 .Random.seed 的两个概率，而不是 0.5：

set.seed(234908)
x <- rbinom(n=1,size=60,prob=0.4)
temp1 <- .Random.seed

set.seed(234908)
x <- rbinom(n=1,size=60,prob=0.3)
temp2 <- .Random.seed

any(temp1!=temp2)
> [1] FALSE

比较 prob=0.5 与 prob!=0.5 的种子向量 .Random.seed：

set.seed(234908)
x <- rbinom(n=1,size=60,prob=0.5)
temp1 <- .Random.seed

set.seed(234908)
x <- rbinom(n=1,size=60,prob=0.3)
temp2 <- .Random.seed
any(temp1!=temp2)
> [1] TRUE

temp1==temp2
> [1]  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE
> [8]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
...

我在所有prob=0.5 与所有其他概率的比较中发现了这一点 在集合 0.1, 0.2, ..., 0.9 中。同样，如果我比较 prob 的任何值 0.1, 0.2, ..., 0.9 除了 0.5，.Random.seed 向量总是逐个元素相等。这些事实也适用于 rbinom() 中的奇数或偶数 size。

为了让它更奇怪（我很抱歉这有点令人费解 - 它与我的函数编写方式有关），当我使用保存为向量中元素的概率时，如果 0.5 是第一个元素，我会遇到同样的问题，但不是第二个。以下是本案例的示例：

第一种情况：0.5是向量中引用的第一个概率

set.seed(234908)
MNAR <- c(0.5,0.3)
x <- rbinom(n=1,size=60,prob=MNAR[1])
y <- rbinom(n=1,size=50,prob=MNAR[2])
temp1 <- .Random.seed

set.seed(234908)
MNAR <- c(0.1,0.3)
x <- rbinom(n=1,size=60,prob=MNAR[1])
y <- rbinom(n=1,size=50,prob=MNAR[2])
temp2 <- .Random.seed

any(temp1!=temp2)
> [1] TRUE

any(temp1!=temp2)
> [1]  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE
> [8]  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE

第二种情况：0.5是向量中引用的第二个概率

set.seed(234908)
MNAR <- c(0.3,0.5)
x <- rbinom(n=1,size=60,prob=MNAR[1])
y <- rbinom(n=1,size=50,prob=MNAR[2])
temp1 <- .Random.seed

set.seed(234908)
MNAR <- c(0.1,0.3)
x <- rbinom(n=1,size=60,prob=MNAR[1])
y <- rbinom(n=1,size=50,prob=MNAR[2])
temp2 <- .Random.seed

any(temp1!=temp2)
> [1] FALSE

我再次发现，尽管prob 和size 使用了这些值，但这种模式仍然成立。谁能给我解释一下这个谜？这造成了很大的问题，因为应该相同的结果会出现不同的结果，因为在 prob=0.5 但在其他情况下，由于某种原因，种子的使用/计算方式不同。

【问题讨论】：

这种行为的一个更简洁的例子可能是set.seed(123);rbinom(1,60,0.5);rbinom(1,60,0.3); set.seed(123);rbinom(1,60,0.2);rbinom(1,60,0.3); set.seed(123);rbinom(1,60,0.4);rbinom(1,60,0.3) ?

我会去svn.r-project.org/R/trunk/src/nmath/rbinom.c，搜索unif_rand()，然后按照逻辑...

使用上面@joran 的例子的一个有趣的结果是，如果你让prob = 0.2 或prob = 0.4 画两个数字而不是一个数字，你就会重新站起来。它表明prob = 0.5 需要绘制的随机数是其他概率的两倍。该理论还通过在 OP 的 x <- rbinom(n=1,size=60,prob=0.3) 案例中将 60 替换为 120 来验证。

哈！查看 Ben 所指的代码：n*p >= 30 和 n*p < 30 有部分。前者使用两次对unif_rand() 的调用，后者使用一次调用。现在请注意，您的示例使用了prob = 0.5 和size = 60，即n*p == 30！使用size = 59 进行测试，行为消失！

@flodel：把它写下来作为答案！

【参考方案1】：

所以让我们把我们的 cmets 变成一个答案。感谢 Ben Bolker 为我们提供了代码链接：https://svn.r-project.org/R/trunk/src/nmath/rbinom.c，并提出了追踪 unif_rand() 被调用位置的建议。

快速扫描一下，代码似乎分为两部分，由 cmets 分隔：

/*-------------------------- np = n*p >= 30 : ------------------- */

和

/*---------------------- np = n*p < 30 : ------------------------- */

在其中的每一个中，对unif_rand 的调用次数都不相同（两个对一个。）

因此，对于给定的 size (n)，您的随机种子最终可能会处于不同的状态，具体取决于 prob (p) 的值：无论是否为 size * prob >= 30。

考虑到这一点，您通过示例获得的所有结果现在应该是有意义的：

# these end up in the same state
rbinom(n=1,size=60,prob=0.4) # => np <  30
rbinom(n=1,size=60,prob=0.3) # => np <  30

# these don't
rbinom(n=1,size=60,prob=0.5) # => np >= 30
rbinom(n=1,size=60,prob=0.3) # => np <  30

# these don't
rbinom(n=1,size=60,prob=0.5)  # np >= 30
 rbinom(n=1,size=50,prob=0.3) # np <  30
rbinom(n=1,size=60,prob=0.1)  # np <  30
 rbinom(n=1,size=50,prob=0.3) # np <  30

# these do
rbinom(n=1,size=60,prob=0.3)  # np <  30
 rbinom(n=1,size=50,prob=0.5) # np <  30
rbinom(n=1,size=60,prob=0.1)  # np <  30
 rbinom(n=1,size=50,prob=0.3) # np <  30

【讨论】：

太好了。但是，我想我仍然有麻烦，因为样本大小是在我的代码中随机抽取的（我只是在这里选择了一些数字进行说明）。至少知道哪里出了问题是件好事，但我的模拟仍然会受到影响。这还是让我不舒服，哈哈。谢谢大家！

@Meg 现在我们知道发生了什么，您应该能够通过确保在比较运行中使用相同的prob 值或通过调用set.seed 来处理它代码中的多个（匹配）位置以“重置”公共序列。

为什么不作为功能请求发布到 R-devel？无论参数值如何，似乎 rbinom 都应该对随机数流产生一致的影响； runif 可以。

np>30 算法使用拒绝方法——该方法的每次迭代抽取两个均匀随机样本，这意味着抽取的样本数量将取决于不仅仅是参数，还有RNG的当前状态。我认为这强化了@BrianDigg 的回答，即您想要的行为可能不切实际。 （@MartinMorgan 评论说runif 对流有一致的影响，因此所有r* 应该是不公平的：runif 是基本构建块并直接调用unif_rand。

[续] 除非您始终可以使用确定性转换方法，否则您将无法做出您想要的保证。

【参考方案2】：

我将在这个问题上采取相反的立场，并声称期望是不合适的，并且没有得到文档的支持。该文档没有说明调用 rbinom 会产生什么副作用（特别是在 .Random.seed 上），或者这些副作用在各种情况下可能相同或不同。

rbinom 具有三个参数：n、size 和 prob。您的期望是，对于调用 rbinom 之前的随机种子集，对于给定的 n 和 any 值 size 和 @ 在调用 rbinom 之后，.Random.seed 将是相同的987654333@（或者可能是size 和prob 的任何有限值）。您当然意识到n 的不同值会有所不同。 rbinom 不保证或暗示这一点。

不知道函数的内部，这是无法知道的；正如other answer 所示，算法基于size 和prob 的乘积而有所不同。而且内部结构可能会发生变化，因此这些具体细节可能会发生变化。

至少，在这种情况下，每次调用 rbinom 后得到的 .Random.seed 将是相同的，其中 相同 n、size 和 em>prob。我可以构建一个病态函数，这甚至不是真的：

seedtweak <- function() 
  if(floor(as.POSIXlt(Sys.time())$sec * 10) %% 2) 
    runif(1)
  
  invisible(NULL)

基本上，这个函数看一个时间的十分之一是奇数还是偶数来决定是否抽取一个随机数。运行这个函数，.Random.seed 可能会也可能不会改变：

rs <- replicate(10, 
  set.seed(123) 
  seedtweak()
  .Random.seed
)
all(apply(rs, 1, function(x) Reduce(`==`, x)))

您可以（应该？）希望的最好结果是，输入/参数相同（包括种子）的给定代码集所有将始终给出相同的结果。当只有大多数（或只有部分）参数相同时期望相同​​的结果是不现实的，除非所有调用的函数都做出这些保证。

【讨论】：

+1。我完全同意。此外，要保证这种再现性可能根本不可能。例如，考虑拒绝抽样。

伟大的论点并写出来！