Отбор проб из резервуара

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

Мотивация

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

Просто: алгоритм R

Простой и популярный, но медленный алгоритм Алгоритм R был создан Джеффри Виттером.

Инициализировать массив $R$ проиндексировано из $1$ к $k$, содержащий первые k элементов входных данных $x_{1},...,x_{k}$. Это резервуар .

Для каждого нового входа $x_{i}$, сгенерируйте случайное число j равномерно по $\{1,...,i\}$. Если $j\in \{1,...,k\}$, затем установите $R[j]:=x_{i}.$ В противном случае отбросить $x_{i}$.

Возвращаться $R$ после обработки всех входных данных.

Этот алгоритм работает методом индукции по $i\geq k$.

Доказательство

Когда $i=k$Алгоритм R возвращает все входные данные, обеспечивая тем самым основу для доказательства методом математической индукции.

Здесь индукционная гипотеза заключается в том, что вероятность того, что конкретный входной сигнал будет включен в резервуар непосредственно перед $(i+1)$-й ввод обрабатывается ${\frac {k}{i}}$ и мы должны показать, что вероятность того, что конкретный входной ресурс включен в резервуар, равна ${\frac {k}{i+1}}$ сразу после $(i+1)$-й ввод обрабатывается.

Примените алгоритм R к $(i+1)$-й вход. Вход $x_{i+1}$ включено с вероятностью ${\frac {k}{i+1}}$ по определению алгоритма. Для любого другого ввода $x_{r}\in \{x_{1},...,x_{i}\}$, по гипотезе индукции, вероятность того, что он окажется в резервуаре непосредственно перед $(i+1)$-й ввод обрабатывается ${\frac {k}{i}}$. Вероятность того, что он все еще будет включен в резервуар после $x_{i+1}$ обрабатывается (другими словами, что $x_{r}$ не заменяется на $x_{i+1}$) является ${\frac {k}{i}}\times \left(1-{\frac {1}{i+1}}\right)={\frac {k}{i+1}}$. Последнее следует из предположения, что целое число $j$ генерируется равномерно случайным образом; как только станет ясно, что замена действительно произойдет, вероятность того, что $x_{r}$ в частности заменяется на $x_{i+1}$ является ${\frac {k}{i+1}}{\frac {1}{k}}={\frac {1}{i+1}}$.

Мы показали, что вероятность того, что новый входной сигнал попадет в резервуар, равна вероятности того, что существующий входной сигнал в резервуаре сохранится. Таким образом, по принципу математической индукции мы приходим к выводу, что алгоритм R действительно создает случайную выборку входных данных.

Хотя этот алгоритм концептуально прост и понятен, он должен генерировать случайное число для каждого элемента входных данных, включая элементы, которые отбрасываются. алгоритма равно асимптотическое время работы Таким образом, $O(n)$. Генерация такого количества случайности и линейного времени выполнения приводит к тому, что алгоритм становится излишне медленным, если входная совокупность велика.

Это алгоритм R, реализованный следующим образом:

(* 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]
  end

  // replace elements with gradually decreasing probability
  for i := k+1 to n
    (* randomInteger(a, b) generates a uniform integer from the inclusive range {a, ..., b} *)
    j := randomInteger(1, i)
    if j <= k
        R[j] := S[i]
    end
  end
end

Оптимально: алгоритм L

Если мы сгенерируем $n$ случайные числа $u_{1},...,u_{n}\sim U[0,1]$ независимо, то индексы наименьшего $k$ из них представляет собой однородную выборку k-подмножеств $\{1,...,n\}$.

Этот процесс можно выполнить, не зная $n$:

Держите самый маленький $k$ из $u_{1},...,u_{i}$ что было замечено до сих пор, а также $w_{i}$, индекс крупнейшего среди них. Для каждого нового $u_{i+1}$, сравните его с $u_{w_{i}}$. Если $u_{i+1}<u_{w_{i}}$, затем отбросить $u_{w_{i}}$, магазин $u_{i+1}$, и установите $w_{i+1}$ быть индексом крупнейшего среди них. В противном случае отбросить $u_{i+1}$, и установите $w_{i+1}=w_{i}$.

Теперь соедините это с потоком входных данных. $x_{1},...,x_{n}$. Каждый раз, когда некоторые $u_{i}$ принято, сохраните соответствующий $x_{i}$. Каждый раз, когда некоторые $u_{i}$ отбрасывается, отбрасываем соответствующий $x_{i}$.

Этот алгоритм все еще нуждается $O(n)$ случайные числа, таким образом принимая $O(n)$ время. Но это можно упростить.

Первое упрощение: нет необходимости тестировать новые $u_{i+1},u_{i+2},...$ один за другим, поскольку вероятность того, что следующая приемка произойдет в $u_{i+l}$ является $(1-u_{w_{i}})^{l-1}\,u_{w_{i}}=(1-u_{w_{i}})^{l-1}-(1-u_{w_{i}})^{l}$, то есть интервал $l$ принятия следует геометрическому распределению .

Второе упрощение: необязательно запоминать весь массив наименьших $k$ из $u_{1},...,u_{i}$ это было замечено до сих пор, но просто $w$, самый крупный среди них. Это основано на трех наблюдениях:

• Каждый раз что-то новое $x_{i+1}$ выбран для ввода в хранилище, равномерно случайная запись в хранилище отбрасывается.
• $u_{i+1}\sim U[0,w]$
• $w_{i+1}$ имеет то же распределение, что и $\max\{u_{1},...,u_{k}\}$, где все $u_{1},...,u_{k}\sim U[0,w]$ независимо. Это может быть отобрано путем первой выборки. $u\sim U[0,1]$, затем берём $w\cdot u^{\frac {1}{k}}$.

Это алгоритм L , который реализуется следующим образом:

(* 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]
  end

  (* random() generates a uniform (0,1) random number *)
  W := exp(log(random())/k)

  while i <= n
      i := i + floor(log(random())/log(1-W)) + 1
      if i <= n
          (* replace a random item of the reservoir with item i *)
          R[randomInteger(1,k)] := S[i]  // random index between 1 and k, inclusive
          W := W * exp(log(random())/k)
      end
  end
end

Этот алгоритм вычисляет три случайных числа для каждого элемента, который становится частью резервуара, и не тратит время на элементы, которые этого не делают. Таким образом, его ожидаемое время работы равно $O(k(1+\log(n/k)))$, что является оптимальным. В то же время его легко реализовать и он не зависит от случайных отклонений от экзотических или трудновычислимых распределений.

Со случайной сортировкой

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

(*
  S is a stream of items to sample
  S.Current returns current item in stream
  S.Next advances stream to next position
  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..?])
  H := new min-priority-queue
  while S has data
    r := random()   // uniformly random between 0 and 1, exclusive
    if H.Count < k
      H.Insert(r, S.Current)
    else
      // keep k items with largest associated keys
      if r > H.Minimum
        H.Extract-Min()
        H.Insert(r, S.Current)
    end
    S.Next
  end
  return items in H
end

Ожидаемое время работы этого алгоритма равно $O(n+k\log k\log(n/k))$ и он актуален главным образом потому, что его можно легко распространить на предметы с весом.

Взвешенная случайная выборка

Этот метод, называемый также последовательной выборкой, некорректен в том смысле, что он не позволяет получить априорно фиксированные вероятности включения. Некоторые приложения требуют, чтобы вероятности выборки элементов соответствовали весам, связанным с каждым элементом. Например, может потребоваться выборка запросов в поисковой системе с весом, определяющим количество раз, когда они выполнялись, чтобы можно было проанализировать выборку на предмет общего влияния на взаимодействие с пользователем. Пусть вес предмета i будет $w_{i}$, а сумма всех весов равна W . Существует два способа интерпретации весов, присвоенных каждому элементу набора:

1. В каждом раунде вероятность того, что каждый невыбранный элемент будет выбран в этом раунде, пропорциональна его весу относительно весов всех невыбранных элементов. Если X — текущая выборка, то вероятность элемента $i\notin X$ быть выбранным в текущем раунде $\textstyle w_{i}/(W-\sum _{j\in X}{w_{j}})$.
2. Вероятность включения каждого предмета в случайную выборку пропорциональна его относительному весу, т. е. $w_{i}/W$. Обратите внимание, что эта интерпретация может быть недостижима в некоторых случаях, например, $k=n$.

Алгоритм A-Res

Следующий алгоритм был предложен Эфраимидисом и Спиракисом и использует интерпретацию 1:

(*
  S is a stream of items to sample
  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..?])
  H := new min-priority-queue
  while S has data
    r := random() ^ (1/S.Weight)   // random() produces a uniformly random number in (0,1)
    if H.Count < k
      H.Insert(r, S.Current)
    else
      // keep k items with largest associated keys
      if r > H.Minimum
        H.Extract-Min()
        H.Insert(r, S.Current)
      end
    end
    S.Next
  end
  return items in H
end

Этот алгоритм идентичен алгоритму, приведенному в разделе «Отбор проб из резервуара со случайной сортировкой», за исключением генерации ключей элементов. Алгоритм эквивалентен присвоению каждому элементу ключа. $r^{1/w_{i}}$ где r — случайное число, а затем выбирают k элементов с наибольшими ключами. Эквивалентно, более численно стабильная формулировка этого алгоритма вычисляет ключи как $-\ln(r)/w_{i}$ и выберите k предметов с наименьшими ключами.

Алгоритм A-ExpJ

Следующий алгоритм является более эффективной версией A-Res , также предложенной Эфраимидисом и Спиракисом:

(*
  S is a stream of items to sample
  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 the priority queue
    Minimum -> minimum key of any item in the priority queue
    Extract-Min() -> Remove the item with minimum key
    Insert(Key, Item) -> Adds item with specified key
 *)
ReservoirSampleWithJumps(S[1..?])
  H := new min-priority-queue
  while S has data and H.Count < k
    r := random() ^ (1/S.Weight)   // random() produces a uniformly random number in (0,1)
    H.Insert(r, S.Current)
    S.Next
  end
  X := log(random()) / log(H.Minimum) // this is the amount of weight that needs to be jumped over
  while S has data
    X := X - S.Weight
    if X <= 0
      t := H.Minimum ^ S.Weight
      r := random(t, 1) ^ (1/S.Weight) // random(x, y) produces a uniformly random number in (x, y)
    
      H.Extract-Min()
      H.Insert(r, S.Current)

      X := log(random()) / log(H.Minimum)
    end
    S.Next
  end
  return items in H
end

Этот алгоритм использует те же математические свойства, которые используются в A-Res , но вместо вычисления ключа для каждого элемента и проверки, следует ли вставить этот элемент или нет, он вычисляет экспоненциальный переход к следующему элементу, который будет вставлен. Это позволяет избежать необходимости создавать случайные варианты для каждого элемента, что может быть дорогостоящим. Количество необходимых случайных величин уменьшено с $O(n)$ к $O(k\log(n/k))$ в ожидании, где $k$ - размер резервуара, и $n$ количество элементов в потоке.

Алгоритм А-Чао

Предупреждение: следующее описание неверно, см. оригинальную статью Чао и обсуждение здесь .

Следующий алгоритм, предложенный М.Т. Чао, использует интерпретацию 2: и Тилле (2006).

(*
  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
  end
  for i := k+1 to n
    WSum := WSum + S[i].Weight
    p := S[i].Weight / WSum // probability for this item
    j := random();          // uniformly random between 0 and 1
    if j <= p               // select item according to probability
        R[randomInteger(1,k)] := S[i]  //uniform selection in reservoir for replacement
    end
  end
end

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

Обратите внимание, что Чао не указывает, как отбирать первые k элементов. Он просто предполагает, что у нас есть какой-то другой способ выбирать их пропорционально их весу. Чао: «Предположим, что у нас есть план выборки фиксированного размера по отношению к S_k в момент времени A первого порядка ; такой, что вероятность включения X_t равна π(k; i) ».

Алгоритм А-Чао с переходами

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

(*
  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
  end
  j := random() // uniformly random between 0 and 1
  pNone := 1 // probability that no item has been selected so far (in this jump)
  for i := k+1 to n
    WSum := WSum + S[i].Weight
    p := S[i].Weight / WSum // probability for this item
    j -= p * pNone
    pNone := pNone * (1 - p)
    if j <= 0
        R[randomInteger(1,k)] := S[i]  //uniform selection in reservoir for replacement
        j = random()
        pNone := 1
    end
  end
end

Связь с перетасовкой Фишера-Йейтса

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

(* S has the input, R will contain the output permutation *)
Shuffle(S[1..n], R[1..n])
  R[1] := S[1]
  for i from 2 to n do
      j := randomInteger(1, i)  // inclusive range
      R[i] := R[j]
      R[j] := S[i]
  end
end

Обратите внимание: хотя остальные карты перетасованы, только первые k в данном контексте важны . Следовательно, массиву R необходимо отслеживать карты только в первых k позициях во время перемешивания, что уменьшает объем необходимой памяти. Усекая R до длины k , алгоритм соответствующим образом модифицируется:

(* S has items to sample, R will contain the result *)
ReservoirSample(S[1..n], R[1..k])
  R[1] := S[1]
  for i from 2 to k do
      j := randomInteger(1, i)  // inclusive range
      R[i] := R[j]
      R[j] := S[i]
  end
  for i from k + 1 to n do
      j := randomInteger(1, i)  // inclusive range
      if (j <= k)
          R[j] := S[i]
      end
  end
end

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

Ограничения

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