高效的递归随机采样

Question

想象一下以下格式的 df：

   ID1 ID2
1    A   1
2    A   2
3    A   3
4    A   4
5    A   5
6    B   1
7    B   2
8    B   3
9    B   4
10   B   5
11   C   1
12   C   2
13   C   3
14   C   4
15   C   5

问题是为 ID1 中的第一个唯一值随机选择一行（最好调整为 n 行），从数据集中删除相应的 ID2 值，从剩余的 ID2 值池中随机选择一个值作为第二个 ID1 值（即递归），等等。

因此，例如，对于第一个 ID1 值，它将执行sample(1:5, 1)，结果为 2。对于第二个 ID1 值，它将执行sample(c(1, 3:5), 1)，结果为 3。对于第三个 ID1 值，它将执行sample(c(1, 4:5), 1)，结果为 5。不可能出现至少没有一个唯一的 ID2 值可分配给特定 ID1 的情况。然而，当选择多个ID2值（例如三个）时，可能会出现数量不够的情况；在这种情况下，请尽可能选择。最后，结果应该具有类似的格式：

  ID1 ID2
1   A   2
2   B   3
3   C   5

它应该足够高效来处理相当大的数据集（ID1 中有数万个唯一值，每个 ID2 中有数十万个唯一值）。

我尝试了多种方法来解决这个问题，但老实说，它们都没有意义，而且可能只会造成混乱，所以我不会在这里分享它们。

样本数据：

df <- data.frame(ID1 = rep(LETTERS[1:3], each = 5),
                 ID2 = rep(1:5, 3))

Answer 1

可能的解决方案

以下是一些方法：

• Reduce使用+进行基本 R 递归subset
• 最大二分匹配使用igraph
• for使用循环的R 基础动态规划

1. 递归

您可以尝试下面的代码（Reduce应用于递归添加未访问的ID2值）

lst <- split(df, ~ID1)
lst[[1]] <- lst[[1]][sample(1:nrow(lst[[1]]), 1), ]
Reduce(
  function(x, y) {
    y <- subset(y, !ID2 %in% x$ID2)
    rbind(x, y[sample(nrow(y), 1), ])
  },
  lst
)

这使

   ID1 ID2
4    A   4
7    B   2
11   C   1

2.二分匹配

我们可以看到，这个问题可以解释为图论中的匹配问题

library(igraph)
library(dplyr)

g <- df %>%
  arrange(sample(n())) %>%
  graph_from_data_frame() %>%
  set_vertex_attr(
    name = "type",
    value = names(V(.)) %in% df$ID1
  )

type.convert(
  setNames(
    rev(
      stack(
        max_bipartite_match(g)$matching[unique(df$ID1)]
      )
    ), names(df)
  ),
  as.is = TRUE
)

我们可以得到

  ID1 ID2
1   A   2
2   B   5
3   C   1

3.for循环动态规划

  lst <- with(df, split(ID2, ID1))
  v <- c()
  for (k in seq_along(lst)) {
    u <- lst[[k]][!lst[[k]] %in% v]
    v <- c(v, u[sample(length(u), 1)])
  }
  type.convert(
    data.frame(ID1 = names(lst), ID2 = v),
    as.is = TRUE
  )

这使

  ID1 ID2
1   A   4
2   B   5
3   C   3

Answer 2

我认为这个算法可以达到你想要的效果，但是效率不是很高。它可以为其他人提供更快解决方案的起点。

all_ID1 <- unique(df$ID1)
available <- unique(df$ID2)
new_ID2 <-  numeric(length(all_ID1))

for(i in seq_along(all_ID1))
{
  ID2_group <- df$ID2[df$ID1 == all_ID1[i]]
  sample_space <- ID2_group[ID2_group %in% available]
  new_ID2[i]<- sample(sample_space, 1)
  available <- available[available != new_ID2[i]]
}

data.frame(ID1 = all_ID1, ID2 = new_ID2)
#>   ID1 ID2
#> 1   A   5
#> 2   B   1
#> 3   C   2

请注意，如果您用完唯一的 ID2 值，这将不起作用。例如，如果 ID1 列中有字母 A:F，每个字母的 ID2 值为 1:5，那么当您为 ID1 值“F”选择 ID2 值时，已不存在唯一的 ID2 值，因为数字 1 到 5 都已分配给字母 A:E。您没有在问题中说明当没有唯一的 ID2 值分配给特定 ID1 时会发生什么 - 它们应该是NA，还是允许重复？

编辑

以下修改允许任意n选择。如果所有可用数字都用完，则样本空间会得到补充：

AC_function <- function(ID1, ID2, n = 1)
{
  all_ID1   <- rep(unique(ID1), each = n)
  available <- unique(ID2)
  new_ID2   <- numeric(length(all_ID1))

   for(i in seq_along(all_ID1))
   {
     ID2_group    <- ID2[ID1 == all_ID1[i]]
     sample_space <- ID2_group[ID2_group %in% available]
     
     if(length(sample_space) < 1) {
        available    <- unique(ID2)
        sample_space <- ID2_group[ID2_group %in% available]
     }
     if(length(sample_space) == 1) {
        new_ID2[i] <- sample_space
        available <- available[available != new_ID2[i]]
     }
     else {
        new_ID2[i]   <- sample(sample_space, 1)
        available    <- available[available != new_ID2[i]]
     }
   }

  data.frame(ID1 = all_ID1, ID2 = new_ID2)
}

例如：

AC_function(df$ID1, df$ID2)
#>   ID1 ID2
#> 1   A   2
#> 2   B   4
#> 3   C   5

AC_function(df$ID1, df$ID2, n = 2)
#>   ID1 ID2
#> 1   A   1
#> 2   A   2
#> 3   B   5
#> 4   B   4
#> 5   C   3
#> 6   C   2

^{由reprex 包于 2021 年 11 月 3 日创建(v2.0.0)}

Answer 3

免责声明：此解决方案假设数据已排列/排序。如果数据没有排序。请先按ID1列排序，然后使用该功能：

还有另一种方法可以在不使用for-loop/ Recursion 甚至更高级别的函数的情况下执行此操作。我们需要注意sampleR 中的函数是向量化的。因此，如果数据框中的所有组的大小相同，或者大小不断增加，那么您可以使用矢量化样本。

n <- 1 # to be sampled from each group
s <- 5 # size of each group - Note that you have to give the minimum size. 
m <- length(unique(df[[1]])) # number of groups.
size <- min(m*n, s) #Total number of sampled data from the dataframe
samples <- sample(s, size)
index <- rep(seq(s), each = n, length = size) * s - s + samples
df[index, ]

这可以写成一个函数：

sub_sample <- function(data, n){
  st <- table(data[[1]])
  s <- min(st)
  m <- length(st) 
  size <- min(m*n, s) 
  samples <- sample(s, size)
  st1 <- rep(c(0, cumsum(head(st,-1))), each = n, length = size)
  index <- st1 + samples
  data[index, ]
}

sub_sample(df, 1)
   ID1 ID2
1    A   1
7    B   2
13   C   3

sub_sample(df, 2)
   ID1 ID2
1    A   1
5    A   5
8    B   3
7    B   2
14   C   4

请注意，当进行子集化时，n=2我们只有 1 组 C 行。为什么？这是因为 C 组有 5 行。但我们已经为 A 组和 B 组使用了 4 个样本。我们只为 C 组保留了 1 个样本。

速度测试：

什么时候n = 1：

Unit: milliseconds
              expr        min         lq      mean     median        uq       max neval
          f_TIC1()  35.682390  41.610310  53.68912  45.618234  49.88343 227.73160   100
          f_TIC2() 151.643959 166.402584 196.51770 179.098992 192.16335 401.36526   100
          f_TIC3()  11.059033  12.268831  14.53906  13.278606  15.38623  23.32695   100
          f_GKi1()  10.725358  11.879908  14.70369  13.108852  17.86946  26.71074   100
          f_GKi2()  10.816891  11.941173  16.55455  12.989787  17.92708 198.44482   100
          f_GKi3()   8.942479   9.950978  14.94726  10.857187  13.35428 171.08643   100
          f_GKi4()   9.085794   9.834465  13.98820  10.666282  13.20658 191.47267   100
 sub_sample(df, 1)   7.878367   8.737534  11.22173   9.508084  14.22219  19.82063   100

当 时n>1，这段代码可以轻松解决这个问题。其他的需要稍微调整，但它们的速度急剧下降。即使在n = group size. 其他大多数人花费的时间太长甚至失败