将集合拆分为 n 个不相等的子集，关键决定因素是子集中的元素聚合并等于预定数量？
Sad*_*own 4 algorithm partitioning r combinatorics
我正在寻找一组数字，并旨在通过集合分区将它们分成子集。如何生成这些子集的决定因素是确保子集中所有元素的总和尽可能接近由预定分布生成的数字。子集的大小不必相同，每个元素只能在一个子集中。我之前已经通过贪心算法（链接在这里）获得了关于这个问题的指导，但我发现集合中的一些较大的数字大大扭曲了结果。因此，我想对这个问题使用某种形式的集合分区。
一个更深层次的潜在问题，我真的很想纠正未来的问题，我发现我被这些类型问题的“蛮力”方法所吸引。（正如您从我下面的代码中看到的，该代码尝试使用折叠通过“蛮力”来解决问题）。这显然是一种完全低效的解决问题的方法，所以我想用更智能的方法来解决这些最小化类型的问题。因此，非常感谢任何建议。
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 天时完全没有答案。碰巧，我认为考虑算法和优化是一个很大的问题，所以感谢发布。
我要注意的第一件事是，您绝对正确，这不是尝试蛮力的问题。您可能会接近正确的答案，但是由于样本和分布点的数量非常多，您将无法找到最佳解决方案。您需要一种迭代方法，仅当元素使拟合更好时才移动元素，并且算法需要在无法使其变得更好时停止。
我的方法是将问题分为三个阶段：
- 切断数据到接近正确仓作为第一近似
- 将元素从有点太大的 bin 移到有点太小的 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. 所以如果我们想要一个函数只返回j给定的分区dist，我们可以这样做：
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 毫秒。