Reservoir sampling

Предположим, мы видим последовательность элементов, но по одному элементу за 1 раз. Мы хотим сохранить один элемент в памяти, и мы хотим, чтобы он был выбран случайным образом из данной последовательности. Если мы знаем общее количество элементов (n), то есть простое решение: выбрать индекс i между 1 и n с равной вероятностью, и выбрать i-ый элемент. Проблема заключается в том, что мы не всегда знаем n заранее. Возможное решение заключается в следующем:

• Сохранить первый элемент в памяти.
• Когда получим ${\displaystyle i}$-й элемент (для ${\displaystyle i>1}$):
  • с вероятностью ${\displaystyle 1/i}$ сохранить новый элемент вместо текущего элемента;
  • с вероятностью ${\displaystyle 1-1/i}$ сохранить текущий элемент и отбросить новый элемент.

Итак:

• когда есть только один элемент, он сохранится с вероятностью 1;
• когда есть 2 элемента, каждый из них сохранится с вероятностью 1/2;
• когда есть 3 элемента, третий элемент сохранится с вероятностью 1/3, а каждый из предыдущих 2 пунктов также сохранится с вероятностью (1/2)(1-1/3) = (1/2)(2/3) = 1/3;
• по индукции легко доказать, что при наличии n элементов, каждый элемент сохранится с вероятностью 1/n.

В своей статье на эту тему наиболее распространенный пример Джеффри Виттер назвал алгоритмом R[1]. Этот простой O(n) алгоритм, как описано в Dictionary of Algorithms and Data Structures[2] состоит из следующих шагов (при условии, что массивы начинаются с индекса 1, а количество элементов, которые нужно выбрать, k, меньше, чем размер исходного массива, S):

/*
  S has items to sample, R will contain the result
*/
ReservoirSample(S[1..n], R[1..k])
  // fill the reservoir array
  for i = 1 to k
      R[i] := S[i]

  // replace elements with gradually decreasing probability
  for i = k+1 to n
    j := random(1, i)   // important: inclusive range
    if j <= k
        R[j] := S[i]

Алгоритм создаёт «резервуар» — массив размера k — и заполняет его первыми k элементами из S. Затем он перебирает все оставшиеся элементы S, пока S не будет исчерпан. На i-ом элементе S алгоритм генерирует случайное число j между 1 и i. Если j меньше или равно k, то j-ый элемент резервуара заменяется на i-й элемент из S. В сущности, для всех i, i-ый элемент S выбирается для включения в резервуар с вероятностью k/i. Аналогичным образом, на каждой итерации j-й элемент резервуара будет заменён с вероятностью 1/k * k/i, то есть, упрощая, с вероятностью 1/i. Можно показать, что, когда алгоритм закончит свою работу, каждый элемент в S имеет одинаковую вероятность (то есть k / длина(S)), что его выберут в резервуар. Чтобы убедиться в этом, рассмотрим следующее доказательство по индукции. После (i-1)-го шага, допустим, вероятность числа попасть в резервуар равна k/(i-1). Поскольку вероятность числа быть заменены на i-ом шаге равна 1/i, то вероятность того, что он выживает на i-ом шаге, равна (i-1)/i. Таким образом, вероятность того, что заданное число попадёт в резервуар после i-го шага, равна произведению этих двух вероятностей, то есть вероятности уже находиться в резервуаре после (i-1)-го шага и вероятности пережить замену на i-ом шаге. Это произведение равно: (k/(i-1)) * ((i-1)/i) = k/i. Следовательно, результат зависит от i, и поэтому утверждение верно по индукции.

Простой алгоритм на основе резервуара может быть разработан с использованием случайной сортировки[3] и реализован с использованием структуры данных приоритетная очередь. Этот алгоритм присваивает случайное число в качестве ключа к каждому элементу и поддерживать k элементов с минимальным значением для ключей. В сущности, это эквивалентно присвоению случайного числа каждому элементу в качестве ключа, последующей сортировке элементов с помощью этих ключей и выборки первых k элементов. В худшем случае время выполнения алгоритма равно ${\displaystyle O(n\log k)}$, а в лучшем случае — ${\displaystyle O(n)}$. Хотя наихудший случай по времени не так хорош, как Алгоритм R, этот алгоритм может быть легко расширен для взвешенной выборки. Заметим, что оба алгоритма могут работать на потоках неизвестной длины.

/*
  S is a stream of items to sample, R will contain the result
  S.Current returns current item in stream
  S.Next advances stream to next position
  max-priority-queue supports:
    Count -> number of items in priority queue
    Maximum() -> returns maximum key value of all items
    Extract-Max() -> Remove the item with max key
    Insert(key, Item) -> Adds item with specified key
*/
ReservoirSample(S[1..?], R[1..k])
  H = new max-priority-queue
  while S has data
    r = Random(0,1)  // important: inclusive range
    if H.Count < k
      H.Insert(r, S.Current)
    else
      if H.Maximum > r
        H.Extract-Max()
        H.Insert(r, S.Current)
    S.Next

Во многих приложениях выборки требуется производить в соответствии с весами, которые присвоены каждому элементу в данном множестве. Например, это может потребоваться для выборки запросов в поисковой машине с весом как количеством раз, когда они были выполнены, так что образец мог бы быть проанализирован на общее воздействие на пользовательский опыт. Есть два способа интерпретировать веса каждого элемента множества:[4]

1. Пусть вес каждого пункта ${\displaystyle w_{i}}$ , а сумма всех весов ${\displaystyle W}$. Мы можем преобразовать вес в вероятность выбора элемента в качестве образца так: ${\displaystyle P_{i}=w_{i}/W}$.
2. Пусть веса двух элементов ${\displaystyle i}$ и ${\displaystyle j}$ будут ${\displaystyle w_{i}}$ и ${\displaystyle w_{j}}$. Пусть вероятность элемент ${\displaystyle i}$ стать выбранным в качестве образца будет ${\displaystyle p_{i}}$, тогда ${\displaystyle p_{j}=p_{i}\max(1,w_{j}/w_{i})}$.

/*
  S is a stream of items to sample, R will contain the result
  S.Current returns current item in stream
  S.Weight  returns weight of current item in stream
  S.Next advances stream to next position
  The power operator is represented by ^
  min-priority-queue supports:
    Count -> number of items in priority queue
    Minimum() -> returns minimum key value of all items
    Extract-Min() -> Remove the item with minimum key
    Insert(key, Item) -> Adds item with specified key
*/
ReservoirSample(S[1..?], R[1..k])
  H = new min-priority-queue
  while S has data
    r = Random(0,1) ^ (1/S.Weight)  // important: inclusive range
    if H.Count < k
      H.Insert(r, S.Current)
    else
      if H.Minimum < r
        H.Extract-Min()
        H.Insert(r, S.Current)
    S.Next

/*
  S has items to sample, R will contain the result
  S[i].Weight contains weight for each item
*/
WeightedReservoir-Chao(S[1..n], R[1..k])
  WSum = 0
  // fill the reservoir array
  for i = 1 to k
      R[i] := S[i]
      WSum = WSum + S[i].Weight
  for i = k+1 to n
    WSum = WSum + S[i].Weight
    p = S[i].Weight / WSum // probability for this item
    j := random(0, 1);     // important: inclusive range
    if j <= p              // select item according to probability
        R[random(1,k)] := S[i]  //uniform selection in reservoir for replacement

Во многих приложениях, объем данных, из которых необходим небольшой образец, слишком велик, и желательно распределить задачи выбора среди многих машин параллельно, чтобы ускорить процесс. Простой подход, который часто используется, хотя он и менее производителен, заключается в том, чтобы присвоить случайное число в качестве ключа каждому элементу, затем выполнить распределенную сортировки и, наконец, получить образец нужного размера из первых k элементов. Если взвешенный образец подходит, тогда ключ вычисляется с использованием ${\displaystyle r^{1/w_{i}}}$ , где ${\displaystyle r}$ — случайное число, а ${\displaystyle w_{i}}$ — вес элемента. Неэффективность этого подхода с очевидностью вытекает из требуемой распределенной сортировки очень больших объемов данных.

Другой, более эффективный подход для распределенной взвешенной случайной выборки заключается в следующем:[6]

1. Распределить данные среди m машин.
2. Каждая машина производит свою собственную взвешенную выборку с использованием ключа ${\displaystyle r^{1/w_{i}}}$ , как описано в предыдущем разделе, и производит выборку размером <= k элементов.
3. Собрать все m образцов размером <= k, получив в итоге общее количество элементов ${\displaystyle n'<=mk}$.
4. Теперь выбрать k элементов из ${\displaystyle n'}$ элементов из шага 3, используя ключ, который уже был вычислен в шаге 2. Это означает, что вместо повторной генерации ключа с использованием генератора случайных чисел в алгоритме выборки, мы используем ключ, который у нас уже вычислен на шаге 2.

Шаг 4 использует ключи из шага 2, потому что мы можем иметь несбалансированное распределение данных по машинам. Например, предположим, что k = 1, машина m1 машина получает только 1 элемент с весом 10, а машина m2 получает 2 элемента, каждый с весом 100. Интуитивно, вероятность элементов из m1 попасть в финальную выборку равна 1/210. В шаге 3, мы получим 1 элемент из m1, как из m2. Если мы пересчитываем ключи в шаге 4, то вероятность того, что элемент из m1 будет в окончательной выборке, составляет 10/110 вместо требуемых 1/210. Теперь заметим, что взвешенная случайная выборка из предыдущего раздела уменьшает максимальное значение ключа в приоритетной очереди по мере обработки всё большего количества элементов. Поэтому, элементы, выбранный на машине с более крупным куском данных, будут иметь низкие значения ключа и, следовательно, более высокий шанс попасть в выборку.

Предположим, кто-то хотел раздать k случайных карт из колоды игральных карт (то есть n=52). Естественным подходом было бы перетасовать колоду и затем взять первые k карт.

В общем случае, перемешивание тоже должно работать, даже если количество карт в колоде заранее не известно, условие которого удовлетворяется вывернутой наизнанку версией тасования Фишер-Йетса:

Чтобы инициализировать массив a из n элементов как случайно перемешанной копии S, индексы обоих массивов начинаются с 0:    a[0] ← S[0]    for i from 1 to n - 1 do        r ← random (0 .. i)        a[i] ← a[r]        a[r] ← S[i]

Обратите внимание, что хотя остальные карты тасуются, только первые k важны в данном контексте.

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

Усекая a до длины k, алгоритм модифицируется соответствующим образом:

Чтобы инициализировать массив a в k случайных элементов из S (длиной n), индексы обоих массивов начинаются с 0:    a[0] ← S[0]    for i from 1 to k - 1 do        r ← random (0 .. i)        a[i] ← a[r]        a[r] ← S[i]

   for i from k to n - 1 do        r ← random (0 .. i)        if (r < k) then a[r] ← S[i]

Так как порядок первых k карт является несущественным, первый цикл может быть удалён, а массив a может быть инициализирован первыми k элементами из S.

Это напоминает алгоритм R.

Ниже приводится простая реализация алгоритма на языке Python, которая выбирает названия страниц из английской Википедии:

import random
SAMPLE_COUNT = 10

# Force the value of the seed so the results are repeatable
random.seed(12345)

sample_titles = []
for index, line in enumerate(open("enwiki-20091103-all-titles-in-ns0")):
        # Generate the reservoir
        if index < SAMPLE_COUNT:
                sample_titles.append(line)
        else:
                # Randomly replace elements in the reservoir
                # with a decreasing probability.
                # Choose an integer between 0 and index (inclusive)
                r = random.randint(0, index)
                if r < SAMPLE_COUNT:
                        sample_titles[r] = line
print sample_titles

Вероятности выбора резервуарных методов обсуждаются в Chao (1982)[5] и Tillé (2006).[7] Если во время первого отбора вероятности равны k/n (или, в случае процедуры Chao, произвольному множеству неравных вероятностей), то во время второго отбора вероятности зависят от порядка, в котором записи отсортированы в оригинальном резервуаре. Проблема преодолевается кубическим методом выборки Deville and Tillé (2004).[8]

Резервуарные выборки делают предположение, что результат помещается в основной памяти, часто подразумевая, что k является постоянной, не зависящей от n. В применениях, где мы хотели бы выбрать большое подмножество входного списка (скажем, треть, то есть k=n/3), нужно использовать другие методы. Были предложены распределенные реализации этой задачи.[9]

• Скользящая средняя
• Тасование Фишер-Йетса

• Very Fast Reservoir Sampling
• Алгоритм Кнута S с примерами на разных языках программирования