Резервуарная выборка

Землетрясения магнитудой 5 и больше случаются около пяти раз в сутки, и в каталоге USGS за 2015–2025 годы их 19 073. Представьте, что они приходят к вам потоком, по одному, и нужно держать честную случайную выборку из десяти: в любой момент каждое уже пришедшее событие должно оказаться в ней с одной и той же вероятностью. Сколько всего придёт, неизвестно, и хранить весь поток нельзя. Нажмите «Пуск»: события ложатся на карту, а резервуар из k штук обновляется по одному простому правилу.

Карта рисует сама себя: землетрясения ложатся вдоль границ литосферных плит. Шанс войти в резервуар у каждого следующего события падает — сначала каждое, потом одно из сотен, к концу одно из двух тысяч, — но поздние события от этого не обделены: их просто намного больше. Кнопка «2000 прогонов» проверяет это прямо: годы в двух тысячах выборок распределены так же, как во всём каталоге, включая 2021-й, когда землетрясений было больше обычного.

Идея

Первые k элементов кладём в резервуар. Элемент номер i > k берём с вероятностью k/i: выбираем случайное j от 1 до i, и если j ≤ k, кладём элемент на место j, вытесняя то, что там было. Иначе элемент пропускаем. В любой момент резервуар — равномерная случайная выборка из всего, что пришло.

Это алгоритм R. Шаг тот же, что в перемешивании Фишера — Йетса, только перемешивание не выписывается целиком: если представить, что каждый новый элемент встаёт на случайное место в перестановке всех пришедших, резервуар — это её первые k мест.

Почему каждый элемент равновероятен

Больше того, равновероятен не только каждый элемент, но и каждый набор из k: все C(n, k) выборок одинаково вероятны. Листинг ниже проверяет это на маленьком примере: при n = 6 и k = 2 все 360 равновероятных цепочек выборов перебираются целиком, и каждый элемент оказывается в резервуаре ровно в 120 из них.

Алгоритм L: прыгать сразу к следующей замене

Алгоритм R тянет случайное число на каждый элемент, хотя поздние почти все пропускаются. Сколько раз резервуар вообще меняется? Элемент номер i что-то заменяет с вероятностью k/i, и ожидание числа замен — сумма этих вероятностей:

$$\sum_{i=k+1}^{n} \frac{k}{i} = k\,(H_n - H_k) \approx k \ln \frac{n}{k}.$$

На потоке землетрясений при k = 10 это 75 замен на 19 073 события. В 1994 году Ким-Хунг Ли предложил алгоритм L, который на пропускаемые элементы не смотрит вовсе: он сразу вычисляет, сколько элементов пропустить до следующей замены. Ключевое наблюдение — порог. Если каждому элементу приписать случайное число из (0, 1) и держать в резервуаре k элементов с наименьшими числами, выборка тоже будет равномерной. Пусть W — наибольшее из этих k чисел. Следующий элемент войдёт в резервуар, только если его число меньше W, то есть каждый следующий элемент — испытание с вероятностью успеха W. Число неудач до первого успеха распределено геометрически, и его можно вытянуть одним случайным числом: ⌊ln u / ln(1 − W)⌋. После замены порог пересчитывается: новый W — это старый, умноженный на u^1/k.

Черты — элементы, которые вошли в резервуар. На логарифмической оси они ложатся равномерно, и это не совпадение: на отрезок от a до b приходится в среднем k·ln(b/a) замен — поровну на каждый десятикратный участок. Миллиард элементов алгоритм L проходит в браузере за доли миллисекунды и тянет несколько сотен случайных чисел вместо миллиарда. Если поток можно перематывать — например, файл с записями одной длины, — пропущенные элементы не нужно даже читать.

Выборка с весами

Иногда элементы неравноправны: частые поисковые запросы важнее редких, сильные землетрясения — слабых. Павлос Эфраимидис и Пол Спиракис в 2006 году предложили для этого простой ключ: элементу с весом w и случайным u из (0, 1) присвоить u^1/w и держать k элементов с наибольшими ключами — в куче, как в любой задаче «k лучших из потока». Первым при этом оказывается элемент номер i с вероятностью ровно w_i/Σw, и так же дальше среди оставшихся. Естественный на вид ключ «u·w» этого не даёт: он предпочитает тяжёлые элементы, но не пропорционально весу. Числа u^1/w при больших весах слипаются у единицы, и сравнивают их по логарифмам; как именно — во второй задаче, на энергиях настоящих землетрясений.

Код

Поток — тот же каталог землетрясений из песочницы. Алгоритм R и алгоритм L с k = 5 на одном генераторе, а в конце — точная проверка равновероятности перебором всех цепочек.

Алгоритм R потратил 19 068 случайных чисел — по одному на каждое событие после первых пяти — и сделал 44 замены при ожидании 40,7. Алгоритм L сделал 45 замен и потратил всего 137 случайных чисел: по три на замену — длина прыжка, место в резервуаре и новый порог — и ещё два, на первый порог и на последний прыжок, который вышел за конец потока. Выборки у них разные, но обе честные.

Сколько это стоит

Алгоритм R — O(n) времени и случайных чисел, O(k) памяти: резервуар и счётчик. Алгоритм L тратит O(k(1 + log(n/k))) случайных чисел и столько же работы, если поток позволяет пропускать элементы, не читая их; если нет, чтение остаётся O(n), но вызовов генератора становится в тысячи раз меньше. Взвешенная выборка — O(n log k) из-за кучи.

Отдельно стоит число замен в резервуаре, k·(H_n − H_k): если резервуар лежит в общей памяти или на диске, каждая замена — запись. За сутки с замером каждую миллисекунду резервуар из 1028 элементов перепишется около 11,7 тысячи раз — один раз на семь с лишним тысяч замеров. Это первая задача в конце статьи.

Где работает

Командная строка. GNU shuf -n k с версии coreutils 8.22 выбирает k случайных строк резервуаром, если вход большой: файл не читается в память целиком. Предложение 2013 года в списке рассылки coreutils так и объясняло выгоду: раз известно, сколько строк нужно на выходе, весь файл держать незачем.

Метрики. Библиотека Dropwizard Metrics хранит замеры времени ответа не все, а выборкой. UniformReservoir — это алгоритм R на 1028 элементах по умолчанию: его квантили описывают всё время жизни программы. ExponentiallyDecayingReservoir — выборка, в которой у свежих замеров вес больше; она построена по методу «прямого затухания» Кормода и соавторов и примерно отражает последние пять минут.

Распределённые выборки. Если поток идёт через несколько серверов, каждый держит свой резервуар и знает, сколько элементов через него прошло. Общая выборка собирается слиянием: сколько элементов взять из каждого резервуара, решает случай с вероятностями, пропорциональными длинам потоков. Тот же приём работает для выборки из нескольких файлов журнала.

История

Задачу о выборке из потока неизвестной длины решали ещё в 1962 году Фан, Мюллер и Резуха. Алгоритм R в нынешнем виде Дональд Кнут привёл во втором томе «Искусства программирования», указав автором Алана Уотермана; независимо его нашли Маклеод и Беллхаус в 1983 году. Джеффри Виттер в 1985 году в статье «Random sampling with a reservoir» дал методу имя и разобрал, как пропускать элементы: его алгоритм Z тратит O(k(1 + log(n/k))) случайных чисел. Ким-Хунг Ли в 1994 году получил ту же оценку гораздо более простым алгоритмом L, а Эфраимидис и Спиракис в 2006-м — взвешенный вариант.

Задачи

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