R 生成的随机数少于使用 rnorm、rexp、rpois 和 runif 指定的数量

Question

我必须为两组大小为 N 的向量生成随机数。

一组的概率是 p，另一组的概率是 q = 1-p。 （例如，对于 p=0.5 的 1000 人口，我必须从一个分布中生成 500 个随机数，从另一个分布中生成 500 个）。 由于这是一个我必须改变'p'的模拟，所以我编写了我的代码来生成这样的：

group1 = rnorm(n = N*p)
group2 = rnorm(n = N*q) # 1st method
group2 = rnorm(n = (N - N*p)) # 2nd method    

使用上述两种方法，R 在 group2 的几行中生成的随机数比它应该生成的随机数少一个（大约 35% 的行使用第一种方法，大约 12% 的行使用第二种方法）。

我也遇到了与 rexp、rpois 和 runif 相同的错误。

以下是两种方法的快照供您参考。

#### EXAMPLE SCRIPT #####

N = 1000
p1 = seq(0.01, 0.99, 0.001)
q1 = 1 - p1


### FIRST METHOD ###

X = data.frame()
for (i in 1:length(p1))
{
X[i, 1] = p1[i]
X[i, 2] = q1[i]
X[i, 3] = length(runif((N * X[i, 1])))
X[i, 4] = length(runif((N * X[i, 2])))
X[i, 5] = X[i, 4] + X[i, 3]
}

table(X[, 5] == 1000) # column three + coulmn four should sum to 1000


### SECOND METHOD ###

Y = data.frame()
for (i in 1:length(p1))
{
Y[i, 1] = p1[i]
Y[i, 2] = q1[i]
Y[i, 3] = length(runif((N * Y[i, 1])))
Y[i, 4] = length(runif((N - N * Y[i, 1])))
Y[i, 5] = Y[i, 3] + Y[i, 4]
}

table(Y[, 5] == 1000) # column three + coulmn four should sum to 1000

Answer 1

R FAQ 7.31 - 舍入错误 - 您的特定问题归结为：

> p=0.32
> p*1000 + (1-p)*1000
[1]1000

嗯，这看起来是正确的。但真的是这样吗？

> (p*1000 + (1-p)*1000) == 1000
[1] FALSE

没有。为什么不？错到什么程度？

> (p*1000 + (1-p)*1000) - 1000
[1] -1.136868e-13

10^-13 中的 1 部分。这意味着：

> length(runif(1000*p))
[1] 320
> length(runif(1000*(1-p)))
[1] 679

因为：

> as.integer(1000*p)
[1] 320
> as.integer(1000*(1-p))
[1] 679

加起来是 999。有关浮点近似的详细信息，请参阅 R FAQ 7.31

解决方法是在处理计数时尽可能使用整数。

> Np = as.integer(1000*p)
> length(runif(Np))
[1] 320
> length(runif(1000-Np))
[1] 680

而不是将q 计算为1-p 并将其乘以N 来尝试得到1000-N*p。