Алгоритм резервуарной выборки

[Akon32]

Шаг N: (для каждого N > K): с вероятностью K/N принимаем решение, выбрать ли данный пришедший элемент. Если решение положительное — перезаписываем в векторе выбранных элементов элемент с индексом N % K (% — остаток от деления) новым значением. Поскольку элемент номер N имеет равные шансы попасть в выборку со всеми предыдущими элементами, то и остаток от деления N % K будет равномерно случайным.

Я правильно понимаю, что если у нас (например) при K=3 последовательность 1,2,3,1,2,3,1,2,3,..., то из-за размещения элементов по правилу N%K мы не имеем шанса получить в выходной выборке больше одной единицы (а также больше одной двойки и одной тройки)?

Как это согласуется с утверждением:

Алгоритм резервуарной выборки позволяет решить эту задачу за O(N) шагов и O(K) памяти. При этом не требуется знать N заранее, а условие случайности выборки ровно K элементов будет чётко соблюдено.

..? Мне кажется, выбор элементов будет частично неслучаен (коллизии элементов на расстоянии K друг от друга), что проявится в некоторых случаях (как в моём примере).

[rafuck]

В примере, похоже, ошибка и 10 (четвертый элемент посл-ти) должна заменить не 13 (первый элемент в резервуаре), а двойку (второй элемент). Метод вычисления индекса на замену (N%K) предполагает, что индексация резервуара начинается с нуля. Но, тем не менее, похоже вы правы, в вашем примере после первых трех шагов получается выборка [1,2,3], а после четвертого шага либо [1,2,3] (новый элемент не попал в выборку), либо [1,1,3] — новый элемент попал в выборку. Вариант [1,2,1], видимо, невозможен.

[mayorovp]

На самом деле алгоритм требует случайного равномерного выбора заменяемого элемента.

[tangro]

Это так. И этот вопрос интересно обсуждается вот в этом видео начиная с момента 41:07. Автор говорит убедительно, но есть один хитрый момент: он утверждает, что random() % N даёт равномерное распределение от 0 до N, а значит нам достаточно его одного и для решения, брать ли элемент и для определения его позиции.

Начать нужно с того, что вообще чёрт его знает, что такое random() — в стандарте С++ есть только rand(), а значит random() — это какой-то их собственный фейсбучный велосипед и не известно, как он работает. Без этого знания было бы действительно надёжнее сгенерировать второе случайное число.

С другой стороны, пример на R в Википедии тоже обходится одним random'oм — но это принципиально другая функция, генерирующая случайное число по равномерному распределению между двумя заданными. Тут действительно, достаточно его одного и для решения брать ли число (J < K) и для использование его в качестве индекса (распределение между 1 и K будет настолько же равномерным, как и между 1 и N).

[mayorovp]

Одного случайного числа и правда достаточно (до тех пор пока N << RAND_MAX), это и без видео понятно.

Вот только правильный алгоритм выглядит так:

1. Взять случайное число X от 1 до N (или от 0 до N-1)
2. Если оно попадает в диапазон от 1 до K (или от 0 до K-1) — заменить выбранный элемент номер X на кандидата номер N

А у вас — вот такой алгоритм получился:

1. Взять случайное число X от 1 до N (или от 0 до N-1)
2. Если оно попадает в диапазон от 1 до K (или от 0 до K-1) — заменить выбранный элемент номер (N%K) на кандидата номер N

Откуда тут взялось число N%K?

[tangro]

Вопрос с N << RAND_MAX решился исправлением кода в первом примере на более корректный — давайте будем считать, что мы генерируем числа в диапазоне от 1 до N именно так :)

На счёт N%K — я хотел упростить, но, видимо, сам себя запутал. Имелось в виду значение X%K, что для рассматриваемого случая X < K вырождается в просто Х.

Сейчас поправлю, спасибо.

[Akon32]

Да, теперь похоже на правду.

[arabesc]

он утверждает, что random() % N даёт равномерное распределение от 0 до N

Это не так, если (RAND_MAX + 1) не кратно N, т.к. числа в диапазоне [0; RAND_MAX%N + 1) будут иметь бОльшую вероятность выпадения, чем [RAND_MAX%N + 1; N).
Например, при RAND_MAX == 2 и N == 2 вероятность выпадения 0 составит 2/3, а для 1 будет 1/3.

p.s. вижу, дальше в комментариях об этом уже написали

[tangro]

если у нас (например) при K=3 последовательность 1,2,3,1,2,3,1,2,3

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

[Akon32]

он гарантирует, что и первая единица, и четвёртая имеют равные шансы попасть в результирующую выборку на любом шаге алгоритма.

Но они не могут попасть в выборку одновременно.
Более наглядный пример: последовательность 1,2,3,4,5,6,7,8,9, K=3.
Невозможно получить выборку, в которой содержатся одновременно элементы 1 и 4 (это конкретный пример, один из многих).
Т.е., в общем виде, при условии наличия в выборке элемента с номером i вероятность наличия в выборке других элементов с номером (i % K) равна 0. Такие элементы, с равными по модулю K индексами, друг друга вытесняют из выборки, но при этом они не вытесняют остальные элементы.

[HappyLynx]

И все же, под взятием K случайных элементов из множества X обычно подразумевается, что если мы возьмем множество всех подмножеств мощности K множества X, то функция вероятности того, что какое-либо из данных подмножеств будет результатом выборки, должна иметь равномерное распределение.
И Akon32 вполне обоснованно показывает, что описанный алгоритм не обеспечивает данной равномерности, т.к. на модельном множестве некоторые из подмножеств элементов (вне зависимости от их значения) мощности K имеют нулевую вероятность оказаться результатом работы алгоритма. Например (при нумерации с 0 и K = 2), вероятность одновременного попадания в результат элементов с номерами 0 и 2 равна нулю.

Цитата из ветки выше:
> Автор говорит убедительно, но есть один хитрый момент: он утверждает, что random() % N даёт равномерное распределение от 0 до N, а значит нам достаточно его одного и для решения, брать ли элемент и для определения его позиции.

Автор очень лихо оставляет незамеченным тот факт, что вид распределения абсолютно ничего не говорит о взаимной независимости случайных величин. А они у него получаются зависимыми, более того, f1(x) = f2(x).

[tangro]

Akon32 был прав и статья была исправлена, чтобы отразить этот момент.

[vibornoff]

Итак, о чём же идёт речь. Выбрать один случайный элемент из вектора — это элементарная задача:

auto result = vect[rand() % vect.size()]; // С++

И это, разумеется, не совсем правда.

Не будем рассматривать тривиальный случай, когда RAND_MAX < vect.size().

Так вот, правдой оно будет тогда и только тогда, когда vect.size() % RAND_MAX == 0.
Для всех остальных случаев значения из диапазона [0 ... vect.size() % RAND_MAX) будут возвращаться немножечко чаще.

Например, в векторе у нас 2 элемента, а RAND_MAX равен 3. Тогда первый элемент вектора будет возвращатся вдвое чаще второго.

[mayorovp]

Наверное, вы имели в виду RAND_MAX % vect.size().

[vibornoff]

Да, именно так. Рад что вы поняли суть.

[arabesc]

Только для равномерного распределения [0; vect.size()) должно выполняться условие (RAND_MAX + 1) % vect.size() == 0, иначе 0 будет выпадать чаще.

[tangro]

Всё так! Спасибо!

Как на счёт вот такого кода:

std::random_device rd;
std::mt19937 gen(rd());
std::uniform_int_distribution<> dis(0, vect.size() — 1);

auto result = vect[dis(gen)];

Это будет работать лучше, верно?

[vibornoff]

Сурово :) Но работать будет.

[tangro]

Заменил код в статье.

[AngReload]

Считается ли такой вариант как K шагов и 2K дополнительной памяти?

JavaScript

function rnd_int_less(n) {
  return Math.floor(Math.random() * n);
}

function reservoir_sampling(k, n_arr) {
  let n = n_arr.length;
  let k_arr = [];
  let j_arr = [];
  // shuffle
  for (let i = 0; i < k; i++) {
    const j = i + rnd_int_less(n - i);
    j_arr.push(j);
    const a = n_arr[i];
    const b = n_arr[j];
    k_arr[i] = b;
    n_arr[i] = b;
    n_arr[j] = a;
  }
  // unshuffle
  for (let i = k - 1; i >= 0; i--) {
    const j = j_arr[i];
    const a = n_arr[i];
    const b = n_arr[j];
    n_arr[i] = b;
    n_arr[j] = a;
  }
  return k_arr;
}