Разбить набор на n неравных подмножеств с ключевым определяющим фактором, состоящим в том, что элементы в подмножестве агрегируют и равны заранее определенной величине?

Question

Я смотрю на набор чисел и собираюсь разбить их на подмножества с помощью разбиения множеств. Решающим фактором того, как будут генерироваться эти подмножества, будет обеспечение того, чтобы сумма всех элементов в подмножестве была как можно ближе к числу, сгенерированному предварительно определенным распределением. Подмножества не должны быть одинакового размера, и каждый элемент может быть только в одном подмножестве. Ранее мне давали руководство по этой проблеме через жадный алгоритм ( Ссылка здесь ), но я обнаружил, что некоторые из больших чисел в наборе резко исказили результаты. Поэтому я хотел бы использовать некоторую форму разбиения множеств для этой проблемы.

Более глубокая основная проблема, которую я действительно хотел бы исправить для будущих проблем, заключается в том, что я нахожу, что меня тянет к подходу "грубой силы" с этим типом проблем. (Как вы можете видеть из моего кода ниже, который пытается использовать сгибы для решения проблемы «грубой силой»). Это, очевидно, совершенно неэффективный способ решения проблемы, и поэтому я хотел бы решить эти проблемы типа минимизации с помощью более интеллектуального подхода в будущем. Поэтому любые советы очень ценятся.

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 будет тем, в каком подмножестве будет каждый элемент (ie, какому размеру он соответствует). Затем я пытался удалить элементы, которые попали в первое подмножество и продолжить оттуда, но я знаю, что мой метод неэффективен.

Заранее спасибо!

Answer 1

Это не тривиальная проблема, так как вы, вероятно, соберетесь из-за полного отсутствия ответов через 10 дней, даже с вознаграждением. Как это бывает, я думаю, что это большая проблема для размышлений об алгоритмах и оптимизациях, поэтому спасибо за публикацию.

Первое, что я хотел бы отметить, это то, что вы абсолютно правы, что это не та проблема с что попробовать грубую силу. Вы можете приблизиться к правильному ответу, но с нетривиальным количеством выборок и точек распространения вы не найдете оптимального решения. Вам нужен итеративный подход, который перемещает элементы, только если они улучшают подбор, а алгоритм должен останавливаться, когда он не может сделать его лучше.

Мой подход здесь состоит в том, чтобы разбить проблему на три этапа :

1. Вырезать данные в приблизительно правильные ячейки в первом приближении
2. Переместить элементы из ячеек, которые тоже немного большой для тех, кто слишком мал. Делайте это итеративно до тех пор, пока больше нет движений, которые оптимизируют ячейки.
3. Поменяйте местами элементы между столбцами для точной настройки подгонки, пока перестановки не станут оптимальными.

Причина для того, чтобы сделать это в таком порядке, состоит в том, что каждый шаг в вычислительном отношении обходится дороже, поэтому вы хотите передать лучшее приближение каждому шагу, прежде чем позволить ему делать свое дело.

Давайте начнем с функции для обрезки данных. примерно в правильные ячейки:

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")

Мы Можно видеть, что подгонка уже довольно хороша при простом разрезе данных, но мы можем сделать намного лучше, переместив элементы соответствующего размера из ячеек большого размера в ящики меньшего размера. Мы делаем это итеративно, пока никакие движения не улучшат нашу посадку. Мы используем два вложенных цикла while, что должно заставить нас беспокоиться о времени вычислений, но мы начали с близкого совпадения, поэтому нам не нужно слишком много ходов до того, как l oop остановится:

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 выборкам. Кажется, достаточно хорошо для большинства приложений. Он будет расти экспоненциально с большими выборками (около 10 секунд на моем P C для 10 000 выборок), но к этому моменту относительная тонкость выборок означает, что cut_elements %>% move_elements уже дает логарифмическую сумму квадратов ниже -30 и следовательно, было бы исключительно хорошее совпадение без точной настройки swap_elements. Это займет около 30 мс с 10 000 выборок.

Answer 2

Чтобы добавить к отличному ответу @AllanCameron, вот решение, которое использует высокоэффективную функцию comboGeneral из RcppAlgos^*.

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 (например, j в случае ОП), такое, что сумма находится в пределах допуска sum(v) * d[i] для некоторого индекса i (d эквивалентно dist в примере ОП). После того, как мы найдем решение a (обратите внимание, что мы ограничиваем количество решений установкой 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