r - R 中 rbinom(prob=0.5) 的不稳定种子行为

Question

我在 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

比较种子向量，.Random.seed ，对于 prob=0.5 与 prob!=0.5:

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向量总是逐个元素相等。这些事实对于奇数或偶数也适用 size rbinom()内.

为了让它变得更奇怪(我很抱歉这有点复杂 - 这与我的函数的编写方式有关)，当我使用保存为向量中的元素的概率时，如果 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 时种子由于某种原因使用/计算不同。但绝无其他情况。

Answer 1

所以让我们将评论转化为答案。感谢 Ben Bolker 通过代码链接让我们走上了正确的道路:https://svn.r-project.org/R/trunk/src/nmath/rbinom.c以及追踪 unif_rand() 调用位置的建议。

快速浏览一下，代码似乎分为两部分，由注释分隔:

/*-------------------------- 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