R生成的随机数比rnorm，rexp，rpois和runif指定的要少 [英] R generates fewer random numbers than specified with rnorm, rexp, rpois and runif

问题描述

我必须为两个大小为N的向量生成随机数。一个组的概率是p，另一个组的概率是q = 1 -p。
（例如，对于p = 0.5的1000人口，我必须从一个分配中产生500个随机数，而从另一个分配500个）。
因为这是一个模拟，我必须改变'p'我写我的代码生成像这样：

pre $ group1 = rnorm（n = N * p）
group2 = rnorm（n = N * q）＃第一种方法
group2 = rnorm（n =（N-N * p））＃第二种方法$使用上述两种方法，R产生的一个随机数少于它在group2的几行中的随机数（b $ b $ / code>

大约有35％的第一行，大约12％的行使用第二种方法）。

我遇到了rexp，rpois和runif的错误

下面是两个方法供您参考的快照。

  ####示例脚本##### 
 
 N = 1000 
 p1 = seq（0.01，0.99，0.001）
 q1 = 1  -  p1 
 
 
 ###第一种方法### 
 
 X = data.frame（）
 for（i in 1：length（p1））
（b）= {b1，b1，b1，b1，b1，b1，b1， X [i，1]）））
 X [i，4] = length（runif（（N * X [i，2]）））
 X [i，5] = X [i，4] + X [i，3] 
} 
 
表（X [，5] == 1000）＃column three + coulmn four sum sum to 1000 
 
 
 ###第二种方法### 
 
 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] 
} 
 
表（Y [，5] == 1000）＃column three + coulmn four sum sum to 1000 
 R常见问题解答7.31  - 舍入错误 - 你特别的问题归结为：  
 
 
 > p = 0.32 
> p 1000 1000 1000 1000 p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p p看起来正确。但真的吗？
 > （p * 1000 +（1-p）* 1000）== 1000 
 [1] FALSE 
  
没有。为什么不？这是怎么回事？
 > （p * 1000 +（1-p）* 1000） -  1000 
 [1] -1.136868e-13 
  
 
 
 <1>在10 ^ -13的部分。这意味着：
 
 
 >长度（runif（1000 * p））
 [1] 320 
> （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）
>长度（runif（Np））
 [1] 320 
> （runif（1000-Np））
 [1] 680 
  
 code> q 作为 1-p ，然后再乘以 N 并得到 1000-N * p 。
 
I have to generate random numbers for two groups of a vector of size N. 

The probability for one group is p, and for the other is q = 1-p. 
(Eg. for a population of 1000 with p=0.5, I have to generate 500 random number from a distribution and 500 from another). 
Since this is a simulation in which I have to vary 'p' I wrote my code to generate like this:
group1 = rnorm(n = N*p)
group2 = rnorm(n = N*q) # 1st method
group2 = rnorm(n = (N - N*p)) # 2nd method    
With both of the above methods, R generates one less random numbers than it should in several rows of group2 (about 35% of rows with the first, and about 12% of rows with the second method).

I run into the same bug with rexp, rpois and runif as well.

Below is the snapshot of both the methods for your reference.
#### 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

解决方案
R FAQ 7.31 - rounding error - your particular problem boils down to this:
> p=0.32
> p*1000 + (1-p)*1000
[1]1000
well that looks correct. But is it really?
> (p*1000 + (1-p)*1000) == 1000
[1] FALSE
No. Why not? How wrong is it?
> (p*1000 + (1-p)*1000) - 1000
[1] -1.136868e-13
1 part in 10^-13. Which means:
> length(runif(1000*p))
[1] 320
> length(runif(1000*(1-p)))
[1] 679
because:
> as.integer(1000*p)
[1] 320
> as.integer(1000*(1-p))
[1] 679
which adds up to 999. See the R FAQ 7.31 for details on floating point approximations

The solution is to work in integers as much as possible when dealing with counts.
> Np = as.integer(1000*p)
> length(runif(Np))
[1] 320
> length(runif(1000-Np))
[1] 680
rather than computing q as 1-p and multiplying that by N to try and get 1000-N*p.

                        
这篇关于R生成的随机数比rnorm，rexp，rpois和runif指定的要少的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持IT屋！