将集合拆分为 n 个不相等的子集，关键决定因素是子集中的元素聚合并等于预定数量？
Split a set into n unequal subsets with the key deciding factor being that the elements in the subset aggregate and equal a predetermined amount?
我正在寻找一组数字，并打算通过集合划分将它们分成子集。如何生成这些子集的决定性因素将确保子集中所有元素的总和尽可能接近预定分布生成的数字。子集的大小不必相同，每个元素只能在一个子集中。我之前曾通过贪婪算法 () 获得有关此问题的指导，但我发现集合中一些较大的数字会严重扭曲结果。因此，我想使用某种形式的集合分区来解决这个问题。
一个更深层次的潜在问题，我真的很想为将来的问题纠正，我发现我被这些类型的问题的“蛮力”方法所吸引。 （正如您从我下面的代码中看到的那样，它尝试使用折叠通过“蛮力”解决问题）。这显然是解决问题的一种完全低效的方法，因此我想用一种更智能的方法来解决这些最小化类型的问题。因此，非常感谢任何建议。
library(groupdata2)
library(dplyr)
set.seed(345)
j <- runif(500,0,10000000)
dist <- c(.3,.2,.1,.05,.065,.185,.1)
s_diff <- 9999999999
for (i in 1:100) {
x <- fold(j, k = length(dist), method = "n_rand")
if (abs(sum(j) * dist[1] - sum(j[which(x$.folds==1)])) < abs(s_diff)) {
s_diff <- abs(sum(j) * dist[1] - sum(j[which(x$.folds==1)]))
x_fin <- x
}
}
这只是一个简化版本，只关注第一个“子集”。 s_diff 将是模拟的理论结果和实际结果之间的最小差异，而 x_fin 将是每个元素所在的子集（即它对应于哪个折叠）。然后我想删除落入第一个子集的元素并从那里继续，但我知道我的方法效率低下。
提前致谢！
这不是一个小问题，因为您可能会在 10 天内从完全没有答案中收集到，即使有赏金。碰巧，我认为这是一个思考算法和优化的大问题，所以感谢您的发帖。
我要指出的第一件事是，您完全正确，这不是那种可以尝试暴力破解的问题。您可能会接近正确答案，但如果样本和分布点数量众多，您将无法找到最佳解决方案。您需要一种迭代方法，仅当它们使拟合更好时才移动元素，并且算法需要在无法使其更好时停止。
我这里的方法是把问题分成三个阶段：
- 将数据剪切到大致正确的 bins 作为第一个近似值
- 将 元素从有点太大的容器移到有点太小的容器中。反复执行此操作，直到不再有任何移动可以优化 bin。
- 交换列之间的元素以微调拟合，直到交换达到最佳。
按此顺序执行的原因是每个步骤的计算成本更高，因此您希望在让每个步骤执行其操作之前为每个步骤传递更好的近似值。
让我们从将数据切割成大致正确的 bin 的函数开始：
cut_elements <- function(j, dist)
{
# Specify the sums that we want to achieve in each partition
partition_sizes <- dist * sum(j)
# The cumulative partition sizes give us our initial cuts
partitions <- cut(cumsum(j), cumsum(c(0, partition_sizes)))
# Name our partitions according to the given distribution
levels(partitions) <- levels(cut(seq(0,1,0.001), cumsum(c(0, dist))))
# Return our partitioned data as a data frame.
data.frame(data = j, group = partitions)
}
我们想要一种方法来评估这个近似值（以及后续的近似值）与我们的答案有多接近。我们可以针对目标分布进行绘图，但使用数字来评估包含在我们的绘图中的拟合优度也很有帮助。在这里，我将使用样本箱和目标箱之间差异的平方和。我们使用日志来使数字更具可比性。数字越小，越合身。
library(dplyr)
library(ggplot2)
library(tidyr)
compare_to_distribution <- function(df, dist, title = "Comparison")
{
df %>%
group_by(group) %>%
summarise(estimate = sum(data)/sum(j)) %>%
mutate(group = factor(cumsum(dist))) %>%
mutate(target = dist) %>%
pivot_longer(cols = c(estimate, target)) ->
plot_info
log_ss <- log(sum((plot_info$value[plot_info$name == "estimate"] -
plot_info$value[plot_info$name == "target"])^2))
ggplot(data = plot_info, aes(x = group, y = value, fill = name)) +
geom_col(position = "dodge") +
labs(title = paste(title, ": log sum of squares =", round(log_ss, 2)))
}
所以现在我们可以做：
cut_elements(j, dist) %>% compare_to_distribution(dist, title = "Cuts only")
我们可以看到，通过简单的数据切割，拟合已经非常好，但是我们可以通过将适当大小的元素从过大的 bin 移动到过小的 bin 来做得更好。我们迭代地这样做，直到没有更多的动作可以提高我们的适应度。我们使用了两个嵌套的 while 循环，这应该让我们担心计算时间，但我们已经开始了一场势均力敌的比赛，所以在循环停止之前我们不应该走太多步：
move_elements <- function(df, dist)
{
ignore_max = length(dist);
while(ignore_max > 0)
{
ignore_min = 1
match_found = FALSE
while(ignore_min < ignore_max)
{
group_diffs <- sort(tapply(df$data, df$group, sum) - dist*sum(df$data))
group_diffs <- group_diffs[ignore_min:ignore_max]
too_big <- which.max(group_diffs)
too_small <- which.min(group_diffs)
swap_size <- (group_diffs[too_big] - group_diffs[too_small])/2
which_big <- which(df$group == names(too_big))
candidate_row <- which_big[which.min(abs(swap_size - df[which_big, 1]))]
if(df$data[candidate_row] < 2 * swap_size)
{
df$group[candidate_row] <- names(too_small)
ignore_max <- length(dist)
match_found <- TRUE
break
}
else
{
ignore_min <- ignore_min + 1
}
}
if (match_found == FALSE) ignore_max <- ignore_max - 1
}
return(df)
}
让我们看看它做了什么：
cut_elements(j, dist) %>%
move_elements(dist) %>%
compare_to_distribution(dist, title = "Cuts and moves")
您现在可以看到匹配非常接近，我们正在努力查看目标数据和分区数据之间是否存在任何差异。这就是为什么我们需要 GOF 的数值测量。
不过，让我们通过在列之间交换 元素以对其进行微调，从而使尽可能 适合。这一步的计算量很大，但我们已经给了它一个近似值，所以它应该没什么可做的：
swap_elements <- function(df, dist)
{
ignore_max = length(dist);
while(ignore_max > 0)
{
ignore_min = 1
match_found = FALSE
while(ignore_min < ignore_max)
{
group_diffs <- sort(tapply(df$data, df$group, sum) - dist*sum(df$data))
too_big <- which.max(group_diffs)
too_small <- which.min(group_diffs)
current_excess <- group_diffs[too_big]
current_defic <- group_diffs[too_small]
current_ss <- current_excess^2 + current_defic^2
all_pairs <- expand.grid(df$data[df$group == names(too_big)],
df$data[df$group == names(too_small)])
all_pairs$diff <- all_pairs[,1] - all_pairs[,2]
all_pairs$resultant_big <- current_excess - all_pairs$diff
all_pairs$resultant_small <- current_defic + all_pairs$diff
all_pairs$sum_sq <- all_pairs$resultant_big^2 + all_pairs$resultant_small^2
improvements <- which(all_pairs$sum_sq < current_ss)
if(length(improvements) > 0)
{
swap_this <- improvements[which.min(all_pairs$sum_sq[improvements])]
r1 <- which(df$data == all_pairs[swap_this, 1] & df$group == names(too_big))[1]
r2 <- which(df$data == all_pairs[swap_this, 2] & df$group == names(too_small))[1]
df$group[r1] <- names(too_small)
df$group[r2] <- names(too_big)
ignore_max <- length(dist)
match_found <- TRUE
break
}
else ignore_min <- ignore_min + 1
}
if (match_found == FALSE) ignore_max <- ignore_max - 1
}
return(df)
}
让我们看看它做了什么：
cut_elements(j, dist) %>%
move_elements(dist) %>%
swap_elements(dist) %>%
compare_to_distribution(dist, title = "Cuts, moves and swaps")
非常接近相同。让我们量化一下：
tapply(df$data, df$group, sum)/sum(j)
# (0,0.3] (0.3,0.5] (0.5,0.6] (0.6,0.65] (0.65,0.715] (0.715,0.9]
# 0.30000025 0.20000011 0.10000014 0.05000010 0.06499946 0.18500025
# (0.9,1]
# 0.09999969
因此，我们有一个非常接近的匹配：每个分区与目标分布的距离小于百万分之一。考虑到我们只有 500 个测量值可以放入 7 个箱子，这真是令人印象深刻。
在检索数据方面，我们没有触及 j 在数据框 df:
中的顺序
all(df$data == j)
# [1] TRUE
并且分区都包含在 df$group 中。因此，如果我们想要一个函数 return 只是给定 dist 的 j 的分区，我们可以这样做：
partition_to_distribution <- function(data, distribution)
{
cut_elements(data, distribution) %>%
move_elements(distribution) %>%
swap_elements(distribution) %>%
`[`(,2)
}
总而言之，我们创建了一种算法，可以创建异常接近的匹配。但是，如果 运行 花费的时间太长，那就不好了。让我们来测试一下：
microbenchmark::microbenchmark(partition_to_distribution(j, dist), times = 100)
# Unit: milliseconds
# expr min lq mean median uq
# partition_to_distribution(j, dist) 47.23613 47.56924 49.95605 47.78841 52.60657
# max neval
# 93.00016 100
仅需 50 毫秒即可拟合 500 个样本。对于大多数应用程序来说似乎已经足够好了。它会随着更大的样本呈指数增长（在我的 PC 上对于 10,000 个样本大约需要 10 秒），但到那时样本的相对精细度意味着 cut_elements %>% move_elements 已经给你一个低于 -30 的对数平方和并且会因此，如果不对 swap_elements 进行微调，这将是一个非常好的匹配。对于 10,000 个样本，这些只需要大约 30 毫秒。
为了补充@AllanCameron 的出色答案，这里有一个解决方案利用了 RcppAlgos*[=59 中的高效函数 comboGeneral =].
library(RcppAlgos)
partDist <- function(v, d, tol_ratio = 0.0001) {
tot_sum <- d * sum(v)
orig_len <- length(v)
tot_len <- d * orig_len
df <- do.call(rbind, lapply(1L:(length(d) - 1L), function(i) {
len <- as.integer(tot_len[i])
vals <- comboGeneral(v, len,
constraintFun = "sum",
comparisonFun = "==",
limitConstraints = tot_sum[i],
tolerance = tol_ratio * tot_sum[i],
upper = 1)
ind <- match(vals, v)
v <<- v[-ind]
data.frame(data = as.vector(vals), group = rep(paste0("g", i), len))
}))
len <- orig_len - nrow(df)
rbind(df, data.frame(data = v,
group = rep(paste0("g", length(d)), len)))
}
我们的想法是找到 v 的一个子集（例如，在 OP 的情况下是 j），使得某些索引 [=] 的总和在 sum(v) * d[i] 的容差范围内21=]（d 相当于 OP 示例中的 dist）。在我们找到 a 解决方案后（N.B。我们通过设置 upper = 1 来限制解决方案的数量），我们将它们分配到一个组中，然后从 v 中删除它们。然后我们进行迭代，直到我们在 v 中只剩下足够的元素，这些元素将被分配给最后一个分配的值（例如 dist[length[dist]].
这是一个使用 OP 数据的示例：
set.seed(345)
j <- runif(500,0,10000000)
dist <- c(.3,.2,.1,.05,.065,.185,.1)
system.time(df_op <- partDist(j, dist, 0.0000001))
user system elapsed
0.019 0.000 0.019
使用@AllanCameron 的绘图函数我们有：
df_op %>% compare_to_distribution(dist, "RcppAlgos OP Ex")
具有相同分布的更大样本呢：
set.seed(123)
j <- runif(10000,0,10000000)
## N.B. Very small ratio
system.time(df_huge <- partDist(j, dist, 0.000000001))
user system elapsed
0.070 0.000 0.071
结果：
df_huge %>% compare_to_distribution(dist, "RcppAlgos Large Ex")
如您所见，解决方案的扩展性非常好。我们可以通过放宽 tol_ratio 以牺牲结果质量为代价来加快执行速度。
对于大型数据集的参考，@AllanCameron 给出的解决方案只用了不到 3 秒，并给出了类似的对数平方和值 (~44)：
system.time(allan_large <- partition_to_distribution(j, dist))
user system elapsed
2.261 0.675 2.938
* 我是 RcppAlgos
的作者
我正在寻找一组数字，并打算通过集合划分将它们分成子集。如何生成这些子集的决定性因素将确保子集中所有元素的总和尽可能接近预定分布生成的数字。子集的大小不必相同，每个元素只能在一个子集中。我之前曾通过贪婪算法 (
一个更深层次的潜在问题，我真的很想为将来的问题纠正，我发现我被这些类型的问题的“蛮力”方法所吸引。 （正如您从我下面的代码中看到的那样，它尝试使用折叠通过“蛮力”解决问题）。这显然是解决问题的一种完全低效的方法，因此我想用一种更智能的方法来解决这些最小化类型的问题。因此，非常感谢任何建议。
library(groupdata2)
library(dplyr)
set.seed(345)
j <- runif(500,0,10000000)
dist <- c(.3,.2,.1,.05,.065,.185,.1)
s_diff <- 9999999999
for (i in 1:100) {
x <- fold(j, k = length(dist), method = "n_rand")
if (abs(sum(j) * dist[1] - sum(j[which(x$.folds==1)])) < abs(s_diff)) {
s_diff <- abs(sum(j) * dist[1] - sum(j[which(x$.folds==1)]))
x_fin <- x
}
}
这只是一个简化版本，只关注第一个“子集”。 s_diff 将是模拟的理论结果和实际结果之间的最小差异，而 x_fin 将是每个元素所在的子集（即它对应于哪个折叠）。然后我想删除落入第一个子集的元素并从那里继续，但我知道我的方法效率低下。
提前致谢！
这不是一个小问题，因为您可能会在 10 天内从完全没有答案中收集到，即使有赏金。碰巧，我认为这是一个思考算法和优化的大问题，所以感谢您的发帖。
我要指出的第一件事是，您完全正确，这不是那种可以尝试暴力破解的问题。您可能会接近正确答案，但如果样本和分布点数量众多，您将无法找到最佳解决方案。您需要一种迭代方法，仅当它们使拟合更好时才移动元素，并且算法需要在无法使其更好时停止。
我这里的方法是把问题分成三个阶段：
- 将数据剪切到大致正确的 bins 作为第一个近似值
- 将 元素从有点太大的容器移到有点太小的容器中。反复执行此操作，直到不再有任何移动可以优化 bin。
- 交换列之间的元素以微调拟合，直到交换达到最佳。
按此顺序执行的原因是每个步骤的计算成本更高，因此您希望在让每个步骤执行其操作之前为每个步骤传递更好的近似值。
让我们从将数据切割成大致正确的 bin 的函数开始：
cut_elements <- function(j, dist)
{
# Specify the sums that we want to achieve in each partition
partition_sizes <- dist * sum(j)
# The cumulative partition sizes give us our initial cuts
partitions <- cut(cumsum(j), cumsum(c(0, partition_sizes)))
# Name our partitions according to the given distribution
levels(partitions) <- levels(cut(seq(0,1,0.001), cumsum(c(0, dist))))
# Return our partitioned data as a data frame.
data.frame(data = j, group = partitions)
}
我们想要一种方法来评估这个近似值（以及后续的近似值）与我们的答案有多接近。我们可以针对目标分布进行绘图，但使用数字来评估包含在我们的绘图中的拟合优度也很有帮助。在这里，我将使用样本箱和目标箱之间差异的平方和。我们使用日志来使数字更具可比性。数字越小，越合身。
library(dplyr)
library(ggplot2)
library(tidyr)
compare_to_distribution <- function(df, dist, title = "Comparison")
{
df %>%
group_by(group) %>%
summarise(estimate = sum(data)/sum(j)) %>%
mutate(group = factor(cumsum(dist))) %>%
mutate(target = dist) %>%
pivot_longer(cols = c(estimate, target)) ->
plot_info
log_ss <- log(sum((plot_info$value[plot_info$name == "estimate"] -
plot_info$value[plot_info$name == "target"])^2))
ggplot(data = plot_info, aes(x = group, y = value, fill = name)) +
geom_col(position = "dodge") +
labs(title = paste(title, ": log sum of squares =", round(log_ss, 2)))
}
所以现在我们可以做：
cut_elements(j, dist) %>% compare_to_distribution(dist, title = "Cuts only")
我们可以看到，通过简单的数据切割，拟合已经非常好，但是我们可以通过将适当大小的元素从过大的 bin 移动到过小的 bin 来做得更好。我们迭代地这样做，直到没有更多的动作可以提高我们的适应度。我们使用了两个嵌套的 while 循环，这应该让我们担心计算时间，但我们已经开始了一场势均力敌的比赛，所以在循环停止之前我们不应该走太多步：
move_elements <- function(df, dist)
{
ignore_max = length(dist);
while(ignore_max > 0)
{
ignore_min = 1
match_found = FALSE
while(ignore_min < ignore_max)
{
group_diffs <- sort(tapply(df$data, df$group, sum) - dist*sum(df$data))
group_diffs <- group_diffs[ignore_min:ignore_max]
too_big <- which.max(group_diffs)
too_small <- which.min(group_diffs)
swap_size <- (group_diffs[too_big] - group_diffs[too_small])/2
which_big <- which(df$group == names(too_big))
candidate_row <- which_big[which.min(abs(swap_size - df[which_big, 1]))]
if(df$data[candidate_row] < 2 * swap_size)
{
df$group[candidate_row] <- names(too_small)
ignore_max <- length(dist)
match_found <- TRUE
break
}
else
{
ignore_min <- ignore_min + 1
}
}
if (match_found == FALSE) ignore_max <- ignore_max - 1
}
return(df)
}
让我们看看它做了什么：
cut_elements(j, dist) %>%
move_elements(dist) %>%
compare_to_distribution(dist, title = "Cuts and moves")
您现在可以看到匹配非常接近，我们正在努力查看目标数据和分区数据之间是否存在任何差异。这就是为什么我们需要 GOF 的数值测量。
不过，让我们通过在列之间交换 元素以对其进行微调，从而使尽可能 适合。这一步的计算量很大，但我们已经给了它一个近似值，所以它应该没什么可做的：
swap_elements <- function(df, dist)
{
ignore_max = length(dist);
while(ignore_max > 0)
{
ignore_min = 1
match_found = FALSE
while(ignore_min < ignore_max)
{
group_diffs <- sort(tapply(df$data, df$group, sum) - dist*sum(df$data))
too_big <- which.max(group_diffs)
too_small <- which.min(group_diffs)
current_excess <- group_diffs[too_big]
current_defic <- group_diffs[too_small]
current_ss <- current_excess^2 + current_defic^2
all_pairs <- expand.grid(df$data[df$group == names(too_big)],
df$data[df$group == names(too_small)])
all_pairs$diff <- all_pairs[,1] - all_pairs[,2]
all_pairs$resultant_big <- current_excess - all_pairs$diff
all_pairs$resultant_small <- current_defic + all_pairs$diff
all_pairs$sum_sq <- all_pairs$resultant_big^2 + all_pairs$resultant_small^2
improvements <- which(all_pairs$sum_sq < current_ss)
if(length(improvements) > 0)
{
swap_this <- improvements[which.min(all_pairs$sum_sq[improvements])]
r1 <- which(df$data == all_pairs[swap_this, 1] & df$group == names(too_big))[1]
r2 <- which(df$data == all_pairs[swap_this, 2] & df$group == names(too_small))[1]
df$group[r1] <- names(too_small)
df$group[r2] <- names(too_big)
ignore_max <- length(dist)
match_found <- TRUE
break
}
else ignore_min <- ignore_min + 1
}
if (match_found == FALSE) ignore_max <- ignore_max - 1
}
return(df)
}
让我们看看它做了什么：
cut_elements(j, dist) %>%
move_elements(dist) %>%
swap_elements(dist) %>%
compare_to_distribution(dist, title = "Cuts, moves and swaps")
非常接近相同。让我们量化一下：
tapply(df$data, df$group, sum)/sum(j)
# (0,0.3] (0.3,0.5] (0.5,0.6] (0.6,0.65] (0.65,0.715] (0.715,0.9]
# 0.30000025 0.20000011 0.10000014 0.05000010 0.06499946 0.18500025
# (0.9,1]
# 0.09999969
因此，我们有一个非常接近的匹配：每个分区与目标分布的距离小于百万分之一。考虑到我们只有 500 个测量值可以放入 7 个箱子，这真是令人印象深刻。
在检索数据方面，我们没有触及 j 在数据框 df:
all(df$data == j)
# [1] TRUE
并且分区都包含在 df$group 中。因此，如果我们想要一个函数 return 只是给定 dist 的 j 的分区，我们可以这样做：
partition_to_distribution <- function(data, distribution)
{
cut_elements(data, distribution) %>%
move_elements(distribution) %>%
swap_elements(distribution) %>%
`[`(,2)
}
总而言之，我们创建了一种算法，可以创建异常接近的匹配。但是，如果 运行 花费的时间太长，那就不好了。让我们来测试一下：
microbenchmark::microbenchmark(partition_to_distribution(j, dist), times = 100)
# Unit: milliseconds
# expr min lq mean median uq
# partition_to_distribution(j, dist) 47.23613 47.56924 49.95605 47.78841 52.60657
# max neval
# 93.00016 100
仅需 50 毫秒即可拟合 500 个样本。对于大多数应用程序来说似乎已经足够好了。它会随着更大的样本呈指数增长（在我的 PC 上对于 10,000 个样本大约需要 10 秒），但到那时样本的相对精细度意味着 cut_elements %>% move_elements 已经给你一个低于 -30 的对数平方和并且会因此，如果不对 swap_elements 进行微调，这将是一个非常好的匹配。对于 10,000 个样本，这些只需要大约 30 毫秒。
为了补充@AllanCameron 的出色答案，这里有一个解决方案利用了 RcppAlgos*[=59 中的高效函数 comboGeneral =].
library(RcppAlgos)
partDist <- function(v, d, tol_ratio = 0.0001) {
tot_sum <- d * sum(v)
orig_len <- length(v)
tot_len <- d * orig_len
df <- do.call(rbind, lapply(1L:(length(d) - 1L), function(i) {
len <- as.integer(tot_len[i])
vals <- comboGeneral(v, len,
constraintFun = "sum",
comparisonFun = "==",
limitConstraints = tot_sum[i],
tolerance = tol_ratio * tot_sum[i],
upper = 1)
ind <- match(vals, v)
v <<- v[-ind]
data.frame(data = as.vector(vals), group = rep(paste0("g", i), len))
}))
len <- orig_len - nrow(df)
rbind(df, data.frame(data = v,
group = rep(paste0("g", length(d)), len)))
}
我们的想法是找到 v 的一个子集（例如，在 OP 的情况下是 j），使得某些索引 [=] 的总和在 sum(v) * d[i] 的容差范围内21=]（d 相当于 OP 示例中的 dist）。在我们找到 a 解决方案后（N.B。我们通过设置 upper = 1 来限制解决方案的数量），我们将它们分配到一个组中，然后从 v 中删除它们。然后我们进行迭代，直到我们在 v 中只剩下足够的元素，这些元素将被分配给最后一个分配的值（例如 dist[length[dist]].
这是一个使用 OP 数据的示例：
set.seed(345)
j <- runif(500,0,10000000)
dist <- c(.3,.2,.1,.05,.065,.185,.1)
system.time(df_op <- partDist(j, dist, 0.0000001))
user system elapsed
0.019 0.000 0.019
使用@AllanCameron 的绘图函数我们有：
df_op %>% compare_to_distribution(dist, "RcppAlgos OP Ex")
具有相同分布的更大样本呢：
set.seed(123)
j <- runif(10000,0,10000000)
## N.B. Very small ratio
system.time(df_huge <- partDist(j, dist, 0.000000001))
user system elapsed
0.070 0.000 0.071
结果：
df_huge %>% compare_to_distribution(dist, "RcppAlgos Large Ex")
如您所见，解决方案的扩展性非常好。我们可以通过放宽 tol_ratio 以牺牲结果质量为代价来加快执行速度。
对于大型数据集的参考，@AllanCameron 给出的解决方案只用了不到 3 秒，并给出了类似的对数平方和值 (~44)：
system.time(allan_large <- partition_to_distribution(j, dist))
user system elapsed
2.261 0.675 2.938
* 我是 RcppAlgos