Резервуарная выборка: k случайных элементов из потока неизвестной длины

Если ты напишешь резервуарную выборку, докажешь, что каждый элемент попадает в неё с вероятностью k/n, и скажешь, что нужно хранить для объединения двух выборок, - листай до «Теперь сам».

Как думаешь, как выбрать одну строку из файла неизвестной длины

Нужно выбрать одну случайную строку из журнала так, чтобы каждая строка имела равные шансы. Журнал читается один раз, строк в нём неизвестно сколько, и хранить можно одну строку.

Ответь себе до чтения: возможно ли это? Кажется, что нет: чтобы выбрать строку с вероятностью 1/n, надо знать n. Возможно. Первую строку берут наверняка, вторую - с вероятностью 1/2 вместо первой, третью - с вероятностью 1/3 вместо той, что сейчас выбрана, и так далее. Когда поток кончится, каждая строка окажется выбранной с вероятностью ровно 1/n.

На что это похоже

Лотерея, в которую участники приходят по одному, и неизвестно, сколько их придёт. В руках у ведущего один приз. Первому он отдаёт приз сразу. Второй подбрасывает монету: выиграл - забирает приз у первого. Третий бросает кубик с тремя гранями, четвёртый - с четырьмя.

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

Механизм: попасть и уцелеть

Резервуар - список из k мест. i - номер текущего элемента, считая с единицы.

для i-го элемента x:
    если i ≤ k:   положить x в резервуар
    иначе:        j = случайное целое от 0 до i − 1
                  если j < k:  резервуар[j] = x

Событие j < k случается с вероятностью k/i: новый элемент попадает в резервуар именно с такой вероятностью и занимает случайное из k мест.

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

Попасть: k/i.

Уцелеть на шаге m, где m > i: новый элемент вытесняет его, если попал в резервуар - вероятность k/m - и занял именно его место из k - вероятность 1/k. Вытеснение случается с вероятностью k/m · 1/k = 1/m, уцелеть - 1 − 1/m = (m − 1)/m.

P(элемент i в выборке) = k/i · i/(i+1) · (i+1)/(i+2) · … · (n−1)/n = k/n

Числитель каждой дроби сокращается со знаменателем предыдущей, и от всего произведения остаются k сверху и n снизу. Номер i исчез: вероятность не зависит от того, когда элемент пришёл. Для первых k элементов то же самое, только «попасть» у них равно единице, а произведение начинается с k/(k+1).

Строго говоря, равные шансы каждого элемента - ещё не всё, что требуется от случайной выборки: нужно, чтобы каждый набор из k элементов был равновероятен. Это тоже выполнено: доказательство - то же перемножение вероятностей, проведённое для набора целиком, и здесь оно опущено.

Правило

Резервуарная выборка размера k: первые k элементов заполняют резервуар; для элемента номер i > k берут j случайно от 0 до i − 1 и при j < k записывают элемент на место j. После n ≥ k элементов каждая позиция потока находится в резервуаре с вероятностью k/n, и все наборы из k позиций равновероятны. При n < k хранят все n элементов; размер равен min(k, n). Память - k мест и счётчик. Чтобы объединить две выборки, нужно знать, сколько элементов видела каждая.

Разбор: один элемент из восьми, шесть тысяч раз

Функция выбирает один элемент потока; random.randrange(videl) == 0 истинно с вероятностью 1/videl. Поток - числа от 0 до 7. Шесть тысяч опытов, и для каждого числа считается, сколько раз выбрали его.

import random

def odin_iz_potoka(potok):
    vybran = -1
    videl = 0
    for x in potok:
        videl = videl + 1
        if random.randrange(videl) == 0:
            vybran = x
    return vybran

random.seed(53)
potok = list(range(8))
skolko = [0] * 8
for opyt in range(6000):
    x = odin_iz_potoka(potok)
    skolko[x] = skolko[x] + 1
print(skolko)

[738, 781, 773, 746, 741, 745, 721, 755]

Ожидается по 750 на каждое число. Все восемь счётчиков лежат между 721 и 781: первый элемент потока, который взят наверняка и должен был пережить семь бросков, выбран 738 раз, последний, у которого один бросок с вероятностью 1/8, - 755 раз.

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

Разбор: три из десяти, и ошибка, которую не видно по одному запуску

Теперь k = 3. Рядом - вариант с правдоподобной ошибкой: новый элемент заменяет случайный элемент резервуара «с вероятностью одна вторая». По полторы тысячи опытов на каждый вариант; печатается, сколько раз каждое число от 0 до 9 оказалось в выборке.

import random

def rezervuar(potok, k):
    vyborka = []
    videl = 0
    for x in potok:
        videl = videl + 1
        if videl <= k:
            vyborka.append(x)
        else:
            j = random.randrange(videl)
            if j < k:
                vyborka[j] = x
    return vyborka

def s_oshibkoj(potok, k):
    vyborka = []
    videl = 0
    for x in potok:
        videl = videl + 1
        if videl <= k:
            vyborka.append(x)
        elif random.random() < 0.5:
            vyborka[random.randrange(k)] = x
    return vyborka

random.seed(54)
potok = list(range(10))
verno = [0] * 10
neverno = [0] * 10
for opyt in range(1500):
    for x in rezervuar(potok, 3):
        verno[x] = verno[x] + 1
    for x in s_oshibkoj(potok, 3):
        neverno[x] = neverno[x] + 1
print(verno)
print(neverno)

[446, 459, 461, 457, 458, 480, 468, 437, 414, 420]
[430, 407, 406, 265, 304, 351, 436, 519, 637, 745]

Правильный вариант: каждое число должно оказаться в выборке 1500 · 3/10 = 450 раз, и счётчики лежат между 414 и 480. Вариант с ошибкой: число 3 попало в выборку 265 раз, число 9 - 745, почти втрое чаще. Постоянная вероятность замены не уравновешивает «попасть» и «уцелеть»: последний элемент потока попадает с вероятностью 1/2 и ни с кем не соревнуется, а четвёртому после своей половины нужно пережить ещё шесть замен.

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

Граница приёма - цена случайных чисел. Алгоритм бросает жребий на каждый элемент, хотя при длинном потоке почти все броски заканчиваются ничем: после миллиона элементов замена происходит один раз на 10⁶/k элементов. Есть вариант, который сразу вычисляет, сколько элементов пропустить до следующей замены; ему нужно порядка k · log(n/k) случайных чисел вместо n. Он полезен, когда элементы можно пропускать не читая - например, перемещаясь по файлу.

Что ломается без этого

Берут первые k элементов. Начало журнала - это ночь, запуск сервиса, прогрев кеша. Выборка из начала потока описывает начало потока.

Берут каждый сотый. Если в данных есть период - запросы чередуются по серверам, записи идут парами, - выборка «каждый N-й» может состоять из записей одного сорта.

Берут каждый элемент с постоянной вероятностью p. Выборка честная, но её размер случаен и растёт вместе с потоком: при неизвестной длине нельзя заранее ни ограничить память, ни гарантировать, что хоть что-то выбрано.

Объединяют выборки без счётчиков. Две выборки по сто строк с серверов, один из которых обработал тысячу строк, а другой миллион, нельзя просто слить и взять сто случайных: строки первого сервера окажутся перепредставлены в тысячу раз.

Зачем это в работе

Отбор примеров из потока, который нельзя ни сохранить, ни перечитать: тысяча случайных запросов за сутки для ручного разбора, выборка событий для отладки, образец строк таблицы для оценки распределения значений. Библиотеки метрик держат резервуар последних измерений времени ответа, чтобы оценивать перцентили по выборке, не храня все измерения.

Командная строка: shuf -n 1000 файл печатает тысячу случайных строк файла (-n - сколько строк взять), не загружая его в память целиком, - внутри резервуар.

Объединение делается так. У каждой выборки есть число виденных элементов, n₁ и n₂. Для выборки размера один: взять элемент первой с вероятностью n₁ / (n₁ + n₂), иначе элемент второй. Для размера k это повторяют k раз: с вероятностью n₁ / (n₁ + n₂) берут случайный ещё не взятый элемент первой выборки и уменьшают n₁ на единицу, иначе - так же из второй, уменьшая n₂. Так распределённая система получает честную выборку из данных, которые ни разу не собирались в одном месте.

Теперь сам

1. Резервуар на один элемент. С какой вероятностью сотый элемент потока заменит текущий выбранный?
2. Резервуар на 10 элементов, прошла тысяча. С какой вероятностью следующий элемент попадёт в резервуар?
3. Резервуар на один элемент, в потоке пять элементов. Перемножь вероятности и найди шанс, что в конце выбран первый.
4. У двух серверов по одной случайной строке; первый видел 300 строк, второй 700. Как выбрать одну случайную строку из всей тысячи?
5. В коде вместо random.randrange(videl) написано random.randrange(k). Что станет с выборкой?