rev*_*er3 5 r sample probability
我有一个关于 R 中样本函数的简单问题。我从长度为 5 的输入向量中随机采样 0 和 1,并将它们相加,该向量指定要运行的试验次数并设置种子以生成可再现的随机数数字。种子按预期工作,但根据我在概率语句中输入的内容,我得到不同的随机数矩阵。在这种情况下,我假设 prob=NULL 应该与 prob=c(0.5,0.5) 相同。为什么不是呢?
vn<-c(12, 44, 9, 17, 28)
> do.call(cbind, lapply(c(1:10),function(X) {set.seed(X); sapply(vn, function(Y) sum(sample(x=c(0,1),size=Y,replace=T)), simplify=TRUE)}))
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 6 7 7 6 6 9 3 6 2 5
[2,] 22 21 20 29 22 24 24 19 25 19
[3,] 4 8 3 5 4 4 4 6 4 2
[4,] 8 4 12 9 11 7 9 10 8 8
[5,] 13 9 11 14 12 14 10 13 11 12
> do.call(cbind, lapply(c(1:10),function(X) {set.seed(X); sapply(vn, function(Y) sum(sample(x=c(0,1),size=Y,replace=T, prob=c(0.5,0.5))), simplify=TRUE)}))
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 6 5 5 6 6 3 9 6 10 7
[2,] 22 23 24 15 22 20 20 25 19 25
[3,] 5 1 6 4 5 5 5 3 5 7
[4,] 9 13 5 8 6 10 8 7 9 9
[5,] 15 19 17 14 16 14 18 15 17 16
Run Code Online (Sandbox Code Playgroud)
更新:
我将采样扩展到 100,并使用输入向量
vn<-seq(0,100,5)
Run Code Online (Sandbox Code Playgroud)
并将没有 prob (test1) 和 prob=c(0.5,0.5) 的输出矩阵的 rowMeans 与预期平均值进行比较。有趣的是,test1 和 test2 的偏差量完全相同,但符号相反。这是为什么?谢谢!
> rowMeans(test1)-seq(0,100,5)/2
[1] 0.00 -0.07 -0.01 -0.35 -0.07 0.19 -0.07 0.24 0.21 0.46 0.20 0.50 -0.37 -0.35 0.00 0.64 -0.59 0.63 -1.19 0.44 -0.38
> rowMeans(test2)-seq(0,100,5)/2
[1] 0.00 0.07 0.01 0.35 0.07 -0.19 0.07 -0.24 -0.21 -0.46 -0.20 -0.50 0.37 0.35 0.00 -0.64 0.59 -0.63 1.19 -0.44 0.38
Run Code Online (Sandbox Code Playgroud)
sample.int正如 Randy 所建议的,根据是否为 NULL使用不同的例程prop。
在您的情况下,它返回相反的结果:
> set.seed(1); sample(c(0,1), size=20, replace=TRUE)
[1] 0 0 1 1 0 1 1 1 1 0 0 0 1 0 1 0 1 1 0 1
> set.seed(1); sample(c(0,1), size=20, replace=TRUE, prob=c(.5,.5))
[1] 1 1 0 0 1 0 0 0 0 1 1 1 0 1 0 1 0 0 1 0
Run Code Online (Sandbox Code Playgroud)
这是怎么回事?
对于前者,我们点击以下行src/main/random.c:546:
for (int i = 0; i < k; i++) iy[i] = (int)(dn * unif_rand() + 1);
Run Code Online (Sandbox Code Playgroud)
这个很简单。 unif_rand()返回 0 到 1 之间的值(并且永远不会返回 1),dn为 2( 中的元素数量x),因此iy[i]设置为1或,具体2取决于是否unif_rand()返回值< .5或>= .5;这是从 中选择的值x。
后者稍微复杂一些。因为prob已指定,do_sample所以调用该函数ProbSampleReplace所以调用处的src/main/random.c:309。revsort这里,概率通过函数at降序排序src/main/sort.c:248。这对概率使用堆排序,并且对于相等概率的二元素向量,它颠倒了顺序。
ProbSampleReplace再次调用,unif_rand()但这次它将其映射到翻转向量顺序后计算的累积概率,因此如果unif_rand()返回一个值,< 0.5则返回第二个值(1在您的示例中)。这是进行映射的代码unif_rand()到 中的值的代码x:
/* compute the sample */
for (i = 0; i < nans; i++) {
rU = unif_rand();
for (j = 0; j < nm1; j++) {
if (rU <= p[j])
break;
}
ans[i] = perm[j];
}
Run Code Online (Sandbox Code Playgroud)
因此,在两个元素的概率相等的情况下,将概率显式设置为c(0.5, 0.5)将返回同一调用的逆函数,而无需设置概率。对于两个以上的元素,它不会总是反转它们,但它不会返回相同的顺序。
这也解释了为什么费尔南多的建议有效。这些值足够接近 0.5,不会更改此示例的结果,并且堆排序按原始顺序返回值。
此表达式返回与第一行代码相同的矩阵:
do.call(cbind, lapply(c(1:10),function(X) {set.seed(X); sapply(vn, function(Y) sum(sample(x=c(1,0),size=Y,replace=T, prob=c(0.5,0.5))), simplify=TRUE)}))
Run Code Online (Sandbox Code Playgroud)
这里,条目的顺序x已反转,以考虑相等值的二元素排序(交换条目)。
当然这都是学术性的。实际上,排列等概率条目的顺序并不重要。
上面的源文件和行号参考R 3.0.2.