这不是一个小问题，因为您可能会在 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 中。因此，如果我们想要一个函数在给定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 个样本。对于大多数应用程序来说似乎已经足够好了。它会随着更大的样本呈指数增长（对于 10,000 个样本，在我的 PC 上大约需要 10 秒），但到那时，样本的相对精细度意味着 cut_elements %>% move_elements 已经为您提供了低于 -30 的平方对数和，因此将是没有swap_elements 微调的非常好的匹配。对于 10,000 个样本，这些只需要大约 30 毫秒。