从数字 1 到 20 生成 100, 000 个大小为 3 的样本，无需放回
Generating 100, 000 samples of size three from numbers 1 to 20, without replacement
我正在尝试从数字 1 到 20 生成大约 100,000 个大小为 3 的样本，无需替换，并在 R 中使用了以下代码：
s <- sample(N,3,pi<-n*x/sum(x),replace=FALSE)
[1] 12 6 17
现在这给了我一个大小为 3 的样本，但是我如何生成其中的 100,000 个？我们还使用了
N<-20 #size of the population we could choose from
n<- 3
x <- runif(N)
pi<-n*x/sum(x)
但我不知道哪里出了问题。任何建议将不胜感激，谢谢。
我们可以通过循环序列
来使用lapply
N1 <- seq_len(100000)
N <- 20
lapply(N1, function(i) sample(N, size =3, replace=FALSE))
你的问题启发了我尝试使用递归采样-with-replacement.
编写一个多采样无替换的实现
让 NS 代表所需样本的数量，NE 代表每个样本的输入集中到 select 的元素数量，我的想法是这可能有益于尽量避免循环 NS sample() 调用，这对于大型 NS 来说会很耗时。相反，我们可以从 运行 单个样本调用开始，采用 NS 值 和 替换，并认为它代表每个样本的 "first selection"。然后，对于每个独特的 selection，我们可以通过 selected 元素减少输入集（和概率加权向量），并递归直到我们达到 NE 级别。通过组合每个（子）样本，我们可以生成一个矩阵，其每一行都将包含来自输入集的 NE 值的无替换样本。
samplesNoReplace <- function(NS,set,NE=length(set),prob=NULL) {
if (NE>1L) {
inds <- sample(seq_along(set),NS,T,prob);
uris <- split(seq_len(NS),inds);
us <- as.integer(names(uris));
res <- base::matrix(set[inds],NS,NE);
for (ui in seq_along(uris)) {
u <- us[ui];
ris <- uris[[ui]];
res[ris,-1L] <- samplesNoReplace(length(ris),set[-u],NE-1L,prob[-u]);
}; ## end for
} else {
res <- base::matrix(sample(set,NS,T,if (length(set)==1L) NULL else prob),ncol=1L);
}; ## end if
res;
}; ## end samplesNoReplace()
演示：
set.seed(10L); samplesNoReplace(10L,1:5,3L,c(10,2,2,2,1));
## [,1] [,2] [,3]
## [1,] 1 3 2
## [2,] 1 4 3
## [3,] 1 2 4
## [4,] 3 2 1
## [5,] 1 3 2
## [6,] 1 4 2
## [7,] 1 4 2
## [8,] 1 2 5
## [9,] 3 1 2
## [10,] 1 2 5
基准测试
library(microbenchmark);
bgoldst <- function() samplesNoReplace(NS,set,NE,prob);
akrun <- function() { N1 <- seq_len(NS); N <- length(set); lapply(N1, function(i) sample(set, size =NE, replace=FALSE,prob)); };
khashaa <- function() { replicate(NS, sample(set, NE,prob=prob), simplify = FALSE); };
## OP's case (100k samples, smallish set, smaller subset)
set.seed(1L);
NS <- 1e5L; set <- 1:20; NE <- 3L; prob <- runif(length(set));
microbenchmark(times=5L,bgoldst(),akrun(),khashaa());
## Unit: milliseconds
## expr min lq mean median uq max neval
## bgoldst() 40.9888 42.69257 46.33044 46.68856 47.40488 53.8774 5
## akrun() 547.3142 564.94249 599.96134 625.07602 631.19658 631.2774 5
## khashaa() 501.1226 521.14871 531.50227 524.65247 549.47600 561.1116 5
## 10k samples, large set, small subset
set.seed(1L);
NS <- 1e4L; set <- 1:1000; NE <- 5L; prob <- runif(length(set));
microbenchmark(times=5L,bgoldst(),akrun(),khashaa());
## Unit: milliseconds
## expr min lq mean median uq max neval
## bgoldst() 2716.1904 2722.8242 2756.9302 2731.2763 2753.5668 2860.7935 5
## akrun() 682.0505 688.3639 691.3169 689.6165 693.9692 702.5842 5
## khashaa() 684.5865 689.2030 698.8313 693.0822 696.1211 731.1638 5
## 1k samples, large set, large subset
set.seed(1L);
NS <- 1e3L; set <- 1:1000; NE <- 500L; prob <- runif(length(set));
microbenchmark(times=1L,bgoldst(),akrun(),khashaa());
## Unit: milliseconds
## expr min lq mean median uq max neval
## bgoldst() 74478.4313 74478.4313 74478.4313 74478.4313 74478.4313 74478.4313 1
## akrun() 350.7270 350.7270 350.7270 350.7270 350.7270 350.7270 1
## khashaa() 353.2574 353.2574 353.2574 353.2574 353.2574 353.2574 1
## 1M samples, small set, necessarily small subset
set.seed(1L);
NS <- 1e6L; set <- 1:4; NE <- 4L; prob <- runif(length(set));
microbenchmark(times=5L,bgoldst(),akrun(),khashaa());
## Unit: milliseconds
## expr min lq mean median uq max neval
## bgoldst() 502.0865 519.1875 602.5631 627.6124 648.3831 715.5459 5
## akrun() 5450.3987 5653.0774 5817.0921 5799.4497 5987.0575 6195.4771 5
## khashaa() 5301.3673 5667.8592 5683.3805 5744.1461 5824.8801 5878.6497 5
## 10M samples, small set, necessarily small subset
set.seed(1L);
NS <- 1e7L; set <- 1:4; NE <- 4L; prob <- runif(length(set));
microbenchmark(times=1L,bgoldst(),akrun(),khashaa());
## Unit: seconds
## expr min lq mean median uq max neval
## bgoldst() 5.023389 5.023389 5.023389 5.023389 5.023389 5.023389 1
## akrun() 75.891354 75.891354 75.891354 75.891354 75.891354 75.891354 1
## khashaa() 69.422056 69.422056 69.422056 69.422056 69.422056 69.422056 1
这个模式很有趣，而且我认为很容易解释。我的函数在许多样本、小集合和小子集上表现出色，因为覆盖所有可能的（子）样本分支所需的递归非常少，而循环解决方案必须迭代并为每个样本调用 sample()。但是我的函数对于更少的样本、大集合和大子集表现严重不佳，因为循环解决方案没有太多迭代要完成，并且（子）样本分支树随着每个新 select 呈指数增长离子。因此，我的函数仅适用于许多样本、小集合和小子集的情况，顺便说一下，这非常准确地描述了您的示例用例。
当然，即使在最不利的时序下，循环解决方案仍然表现不错，大约在我的函数的一个数量级内。此外，在任何情况下都不太可能需要一个小集合的一个小子集的数百万个样本。因此，为了简单起见，我不认为完全忽略此解决方案并始终使用循环方法是不合理的。
我已经尝试了 replicate 命令和 1apply，都为我提供了 1 到 20 的 100,000 个大小为 3 的样本，这很好，但现在我希望能够计算每个数字的频率出现。我知道 9，例如，可能出现 100,000 次，在所有 100,000 个 3-样本中，但更有可能的是，它可能出现大约二十分之一的时间。所以如果我每次有 100,000 个 3 位数的样本，所有数字的总数应该是 300,000，因为为了论证 R 给了我 100,000 个九，每个样本中恰好有 9，那么还剩下二十万个地方对于所有其他数字。我将函数称为 s，并尝试
count1 <- length(which(s == 2)); count1 ，但这说
Error in which(s == 1) : (list) object cannot be coerced to type 'double',
但我不明白那是什么意思。我如何让 R 给我一个所有 1、所有 2 等的准确计数，我假设它们的总数应该是 300,000，因为我们最终在 运行 中得到 300,000 个数字。谢谢。克里斯·莉莉。
我正在尝试从数字 1 到 20 生成大约 100,000 个大小为 3 的样本，无需替换，并在 R 中使用了以下代码：
s <- sample(N,3,pi<-n*x/sum(x),replace=FALSE)
[1] 12 6 17
现在这给了我一个大小为 3 的样本，但是我如何生成其中的 100,000 个？我们还使用了
N<-20 #size of the population we could choose from
n<- 3
x <- runif(N)
pi<-n*x/sum(x)
但我不知道哪里出了问题。任何建议将不胜感激，谢谢。
我们可以通过循环序列
来使用lapply
N1 <- seq_len(100000)
N <- 20
lapply(N1, function(i) sample(N, size =3, replace=FALSE))
你的问题启发了我尝试使用递归采样-with-replacement.
编写一个多采样无替换的实现让 NS 代表所需样本的数量，NE 代表每个样本的输入集中到 select 的元素数量，我的想法是这可能有益于尽量避免循环 NS sample() 调用，这对于大型 NS 来说会很耗时。相反，我们可以从 运行 单个样本调用开始，采用 NS 值 和 替换，并认为它代表每个样本的 "first selection"。然后，对于每个独特的 selection，我们可以通过 selected 元素减少输入集（和概率加权向量），并递归直到我们达到 NE 级别。通过组合每个（子）样本，我们可以生成一个矩阵，其每一行都将包含来自输入集的 NE 值的无替换样本。
samplesNoReplace <- function(NS,set,NE=length(set),prob=NULL) {
if (NE>1L) {
inds <- sample(seq_along(set),NS,T,prob);
uris <- split(seq_len(NS),inds);
us <- as.integer(names(uris));
res <- base::matrix(set[inds],NS,NE);
for (ui in seq_along(uris)) {
u <- us[ui];
ris <- uris[[ui]];
res[ris,-1L] <- samplesNoReplace(length(ris),set[-u],NE-1L,prob[-u]);
}; ## end for
} else {
res <- base::matrix(sample(set,NS,T,if (length(set)==1L) NULL else prob),ncol=1L);
}; ## end if
res;
}; ## end samplesNoReplace()
演示：
set.seed(10L); samplesNoReplace(10L,1:5,3L,c(10,2,2,2,1));
## [,1] [,2] [,3]
## [1,] 1 3 2
## [2,] 1 4 3
## [3,] 1 2 4
## [4,] 3 2 1
## [5,] 1 3 2
## [6,] 1 4 2
## [7,] 1 4 2
## [8,] 1 2 5
## [9,] 3 1 2
## [10,] 1 2 5
基准测试
library(microbenchmark);
bgoldst <- function() samplesNoReplace(NS,set,NE,prob);
akrun <- function() { N1 <- seq_len(NS); N <- length(set); lapply(N1, function(i) sample(set, size =NE, replace=FALSE,prob)); };
khashaa <- function() { replicate(NS, sample(set, NE,prob=prob), simplify = FALSE); };
## OP's case (100k samples, smallish set, smaller subset)
set.seed(1L);
NS <- 1e5L; set <- 1:20; NE <- 3L; prob <- runif(length(set));
microbenchmark(times=5L,bgoldst(),akrun(),khashaa());
## Unit: milliseconds
## expr min lq mean median uq max neval
## bgoldst() 40.9888 42.69257 46.33044 46.68856 47.40488 53.8774 5
## akrun() 547.3142 564.94249 599.96134 625.07602 631.19658 631.2774 5
## khashaa() 501.1226 521.14871 531.50227 524.65247 549.47600 561.1116 5
## 10k samples, large set, small subset
set.seed(1L);
NS <- 1e4L; set <- 1:1000; NE <- 5L; prob <- runif(length(set));
microbenchmark(times=5L,bgoldst(),akrun(),khashaa());
## Unit: milliseconds
## expr min lq mean median uq max neval
## bgoldst() 2716.1904 2722.8242 2756.9302 2731.2763 2753.5668 2860.7935 5
## akrun() 682.0505 688.3639 691.3169 689.6165 693.9692 702.5842 5
## khashaa() 684.5865 689.2030 698.8313 693.0822 696.1211 731.1638 5
## 1k samples, large set, large subset
set.seed(1L);
NS <- 1e3L; set <- 1:1000; NE <- 500L; prob <- runif(length(set));
microbenchmark(times=1L,bgoldst(),akrun(),khashaa());
## Unit: milliseconds
## expr min lq mean median uq max neval
## bgoldst() 74478.4313 74478.4313 74478.4313 74478.4313 74478.4313 74478.4313 1
## akrun() 350.7270 350.7270 350.7270 350.7270 350.7270 350.7270 1
## khashaa() 353.2574 353.2574 353.2574 353.2574 353.2574 353.2574 1
## 1M samples, small set, necessarily small subset
set.seed(1L);
NS <- 1e6L; set <- 1:4; NE <- 4L; prob <- runif(length(set));
microbenchmark(times=5L,bgoldst(),akrun(),khashaa());
## Unit: milliseconds
## expr min lq mean median uq max neval
## bgoldst() 502.0865 519.1875 602.5631 627.6124 648.3831 715.5459 5
## akrun() 5450.3987 5653.0774 5817.0921 5799.4497 5987.0575 6195.4771 5
## khashaa() 5301.3673 5667.8592 5683.3805 5744.1461 5824.8801 5878.6497 5
## 10M samples, small set, necessarily small subset
set.seed(1L);
NS <- 1e7L; set <- 1:4; NE <- 4L; prob <- runif(length(set));
microbenchmark(times=1L,bgoldst(),akrun(),khashaa());
## Unit: seconds
## expr min lq mean median uq max neval
## bgoldst() 5.023389 5.023389 5.023389 5.023389 5.023389 5.023389 1
## akrun() 75.891354 75.891354 75.891354 75.891354 75.891354 75.891354 1
## khashaa() 69.422056 69.422056 69.422056 69.422056 69.422056 69.422056 1
这个模式很有趣，而且我认为很容易解释。我的函数在许多样本、小集合和小子集上表现出色，因为覆盖所有可能的（子）样本分支所需的递归非常少，而循环解决方案必须迭代并为每个样本调用 sample()。但是我的函数对于更少的样本、大集合和大子集表现严重不佳，因为循环解决方案没有太多迭代要完成，并且（子）样本分支树随着每个新 select 呈指数增长离子。因此，我的函数仅适用于许多样本、小集合和小子集的情况，顺便说一下，这非常准确地描述了您的示例用例。
当然，即使在最不利的时序下，循环解决方案仍然表现不错，大约在我的函数的一个数量级内。此外，在任何情况下都不太可能需要一个小集合的一个小子集的数百万个样本。因此，为了简单起见，我不认为完全忽略此解决方案并始终使用循环方法是不合理的。
我已经尝试了 replicate 命令和 1apply，都为我提供了 1 到 20 的 100,000 个大小为 3 的样本，这很好，但现在我希望能够计算每个数字的频率出现。我知道 9，例如，可能出现 100,000 次，在所有 100,000 个 3-样本中，但更有可能的是，它可能出现大约二十分之一的时间。所以如果我每次有 100,000 个 3 位数的样本，所有数字的总数应该是 300,000，因为为了论证 R 给了我 100,000 个九，每个样本中恰好有 9，那么还剩下二十万个地方对于所有其他数字。我将函数称为 s，并尝试 count1 <- length(which(s == 2)); count1 ，但这说 Error in which(s == 1) : (list) object cannot be coerced to type 'double', 但我不明白那是什么意思。我如何让 R 给我一个所有 1、所有 2 等的准确计数，我假设它们的总数应该是 300,000，因为我们最终在 运行 中得到 300,000 个数字。谢谢。克里斯·莉莉。