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



Для таких целей подойдет оценка медианы. То есть такая статистика, что половина значений выборки меньше, а половина больше. Более формально: упорядочим значения выборки X=(x_1,..., x_n) по порядку (x_{[1]}, ..., x_{[n]}) и выберем среди них с порядковым номером floor(n/2). У такой оценки есть несколько преимуществ. Она менее подвержена влиянию ошибочных данных, значение всегда будет из того множества, что встречалось в выборке, но есть и неприятные недостатки, главный из них, это сложность подсчета, даже для довольно распространенных распределений не существует общей формулы расчета (точнее есть, но ее сложно применить на практике, смотрите Распределение порядковой статистики).




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


Но что же делать, если мы не можем позволить себе хранить сто тысяч миллионов значений. С одной стороны можно заняться задачей как отсортировать миллион int-ов в 2 мб оперативки. С другой, в упомянутой выше статье приводится простое решение, которое при некоторых допущениях с некоторой вероятностью приводит к правильному решению. А именно предлагается следующий алгоритм.

Метод Манро-Патерсона


Пусть имеется поток данных длины N (в этом алгоритме, хорошо бы длину потока знать заранее), который мы можем последовательно зачитывать по одному значению. Мы располагаем S ячейками памяти, при этом S<<N и входной поток устроен таким образом, что любая перестановка из N элементов равновероятна. В таком случае, если нам повезет, то алгоритм выдаст значение медианы.

Алгоритм очень прост: хранятся два счетчика (l,h) и подвыборку, состоящую из не более, чем S-2 элементов из входных данных. Счетчик l хранит число элементов, которые меньше минимального в подвыборке, счетчик h — число элементов, которые больше максимального. Сама подвыборка содержит отсортированный список, так что поиск минимума и максимума в нем работает за константу, а вставка элемента линейна по S. Сбалансированные деревья дадут логарифмы для всех операций за небольшую плату дополнительного роста памяти. Когда место под список закончится крайние элементы вытесняются так чтобы уравновесить счетчики l,h (если h>l, то вытесняется минимальный элемент и наоборот).

Так, что после тысячи случайных чисел при S=32, мы можем наблюдать следующую картину.



Если внимательно присмотреться к реализации, становится понятен смысл отсылки к “повезет, не повезет”. Медиана должна попасть в интервал значений хранимых в памяти. В статье приводятся оценки, что это будет происходить довольно часто, если S пропорциональна \sqrt{N}, что в свою очередь означает, что у нас заведомо должна быть информация о длине потока данных и даже в этом случае, мы можем оказаться очень далеко от истинного значения.
Данный алгоритм хорош как первое приближение, так как не очень понятно что же делать, если подсчет не удался и еще не понятно как его распараллеливать.

Дальнейшие исследования пошли в сторону приближенного вычисления медианы и идея приближения очень простая: k-ая порядковая статистика это почти k-1 порядковая статистика или k+1 или какую мы ещё помним недалеко от требуемой. Ошибку в оценке медианы измеряют в доле числа индексов, на которые мы можем “промахиваться”, к n — количеству просмотренных входных данных. Такую ошибку называют \varepsilon n, где \varepsilon — фиксированная допустимая ошибка: например, \varepsilon=0.01 при n=1000 означает, что если бы сохранили всю выборку, то полученная оценка медианы лежала бы между 495-ым и 505-ым отсортированными значениями. В случае если у нас всего 9 значений, то картинка группировки может быть следующей.



Приближенный метод Канна-Гринвальда


Данный метод заимствует идею использования S ячеек памяти, но использует отведенную память для равномерной оценки всех квантилей сразу с заданной \varepsilon n точностью, которая зависит от объема допустимой памяти.

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

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

Для такого дерева можно определить следующие важные операции:
  • Добавление нового значения за log(S);
  • Вычисления медианы и, как легко заметить, произвольного квантиля за линейное по S время;
  • Объединение деревьев при одинаковом опять же за линейное по S время.

В статье приводится (без доказательства) теорема, о том что на случайном наборе данных количество узлов дерева будет расти пропорционально log(S)/\varepsilon, что означает, что алгоритм должен быть очень компактен и спор.

Примеры


В качестве небольшого примера рассмотрим какое-нибудь симметричное и не симметричное распределения с известной медианой. Чтобы далеко не ходить, нормальное и лог-нормальное распределения нам должны подойти. Рассмотрим следующие оценки медианы для выборки в миллион значений:
  • Упорядочивание значений
  • Методом Манро-Патерсона, с параметром равным s=\sqrt{10^7}
  • Методом Канна-Гринвальда c параметром \varepsilon=0.001

Все оценки для каждого из распределений делаются на одних наборах данных и повторяются 100 раз.



В случае нормального распределения среднее и медиана совпадают, в случае лог-нормального распределения они достаточно сильно различаются и использовать одно для оценки другого не имеет смысла (А лучше, так никогда не делать). По парным картинкам видно, что метод полной сортировки очень часто дает совпадение результатов с методом Манро-Патерсона и это правильно, но все-таки есть достаточно много около 8%, как и утверждается в статье, случаев, когда результат Манро-Патерсона отличается от истинного. Метод Канна-Гринвальда дает не плохой, но все же приближенный результат. Ниже видно, что разброс всех методов от истинного значения приблизительно одинаков.



Метод Канна-Гринвальда должен быть достаточно экономичным как по памяти так и по скорости работы, а именно в дереве храниться число элементов пропорциональное логарифму длины входной последовательности и обратно пропорциональное ошибке \varepsilon, а время заполнения пропорционально logN \cdot loglogN.

ПС


Примеры реализации можно найти по ссылке на bitbucket, но это не самая оптимальная реализация.

Уже подготавливая текст я обнаружил, что именно эти статьи и методы упоминаются в курсе алгоритмы обработки потоковых данных в ПОМИ.

Спасибо parpalak за редактор.

Комментарии (6)


  1. SeptiM
    18.08.2015 21:48
    +1

    Интересно, что Канн-Гринвальд детерминированный и всегда дает ответ с нужной погрешностью. Но он мне в свое время показался очень сложным, особенно в смысле понимания, почему работает.

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

    Гарантия доказывается довольно просто через неравенство Чернова.

    P.S. Предлагаю на подобные вещи ставить тэг data streaming.


    1. kokorins
      18.08.2015 22:39

      Надо было упомянуть сэмплирование, но очень хотелось описать алгоритм в котором N не известно заранее. В подходе с сэмплированием описание алгоритма перехода от N к N+1 c объяснением того, что «обновленная» выборка обладает всеми необходимыми свойствами занимает больше места.

      Тег добавил.


  1. vanxant
    18.08.2015 22:39
    +1

    На практике часто случается, что наши данные есть поток каких-то показаний с датчика, причем разрядность датчика часто довольно небольшая (8, 16, 24, пусть даже 32 бита). В этом случае можно посчитать медиану поточно за счет памяти. Просто создаем массив счетчиков, сколько раз встретился 0, сколько 1 и т.д. до maxval, тупо counter[sample]++. По сути строим гистограмму, найти по которой медиану вообще не проблема. При этом гистограмма сама по себе очень полезна в плане статистики, посмотреть распределение, посчитать отклонения, всякие колмогоровы-смирновы, вот это все.


    1. kokorins
      18.08.2015 22:43
      +1

      В тексте есть ссылка на этот подход. И да, практически все пакеты предлагают именно гистограммный подход, например (раз уж я тут пишу на скале), он предлагается в mllib spark-а и breeze (scala backend) для spark. Тут я описываю метод, когда мы хотим обойтись минимальнымы предположениями о структуре входных данных, то есть мы не знаем min/max, не знаем разброс, да и вообще почти ничего не знаем. =)


  1. dedZahar
    19.08.2015 12:00
    +1

    Надо среднюю зарплату так считать


    1. kokorins
      19.08.2015 12:13

      Как раз с этого заметка и возникла. На одном радио был опрос сколько человек получают больше среднего, среди позвонивших было всего 15% таковых. Надеюсь из заметки стало понятнее почему постановка опроса не корректная.