Эффективная оценка медианы

    Итак, у Вас есть какой-то поток данных. Большой такой поток. Или уже готовый набор. И хочется определить какие-то его характеристики. Алгоритм определения минимального и максимального значения могут придумать даже не программисты. Вычисление среднего уже чуть сложнее, но тоже не представляет никаких трудностей — знай подсчитывай себе сумму да инкрементируй счетчик на каждое новое значение. Среднеквадратичное отклонение — все то же самое, только числа другие. А как насчет медианы?

    Для тех, кто забыл, что это такое, напоминаю — медиана (50-й перцентиль) выборки данных — это такое значение, которое делит эту выборку пополам — данные из одной половины имеют значение не меньше медианы, а из второй — не больше. Ценность её заключается в том, что её значение не зависит от величины случайных всплесков, которые могут очень сильно повлиять на среднее.

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

    Сразу предупрежу — предложенный метод обладает рядом ограничений. В частности, он очень плохо работает на отсортированных выборках (но зато очень хорошо работает на более-менее равномерно распределенных). Дальше рассматривается более простой случай неотрицательных значений, для общего случая нужны дополнительные вычисления.

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

    Можно использовать разного рода окна вычисления медианы, например, посчитать точную медиану последних 100 значений, и усреднить с постоянно уменьшающимся весом с предыдущей оценкой. Такой метод будет хорошо работать, но есть гораздо более легкий и практически такой же точный. А именно — просто сдвигать текущую оценку медианы к новому значению на какую-то небольшую константную дельту. В случае, если предыдущая оценка была не точной, то при обработке нового объема данных она приблизится к действительному значению, т.е. станет более точной. А если оценка уже и так точная, то на обработке нового объема данных на половине значений будет сдвиг в одну сторону, а на другой половине — в другую, в итоге оценка вернется к точному значению.

    Важно, что дельта должна иметь одинаковое абсолютное значение для сдвигов как в большую, так и в меньшую сторону, иначе мы получим не 50-й перцентиль. Но теперь встает проблема подбора значения дельты — слишком маленькое даст медленную сходимость, а слишком большое — получим большую погрешность, особенно если дельта сравнима с самим значением медианы. Автоматическое вычисление дельты уже лишает её звания константы, но это и неважно, если в итоге мы получим лучший результат.

    Имеет смысл делать дельту постоянно уменьшающейся, чтобы улучшить сходимость. Но не слишком быстро, иначе, при неблагоприятных условиях оценка рискует никогда не догнать действительное значение медианы. Коэфициент 1/count подходит идеально — он легко вычисляется, и ряд sum(1/n) — расходящийся, так что в конце-концов оценка достигнет действительной медианы, какой бы плохой она ни была изначально.

    Привязывать значение дельты к текущей оценке медианы — рискованно. В неудачных условиях слишком заниженной оценке будет сложно расти. Лучше всего вычислять дельту исходя из другой вполне устойчивой характеристики выборки — среднего значения. С коэфициентом 1/count значение дельты будет постоянно уменьшаться (за исключением немногочисленных случаев, когда текущее среднее вырастет слишком сильно по сравнению с предыдущим). Для данных, которые могут принимать не только неотрицательные значение, такой подход уже не сработает, но мы можем легко выйти из положения, если считать два средних — положительных и отрицательных значений, и использовать сумму их абсолютных значений.

    Итоговый алгоритм оценки медианы почти такой же простой, как и алгоритм вычисления среднего значения:
    sum += value
    count ++
    delta = sum / count / count        // delta = average/count
    if value < median {
      median -= delta
    } else {
      median += delta
    }
    


    Погрешность и скорость сходимости, по моему мнению, у него отличная — на выборке в 100 значений погрешность обычно около 10-20%, на 1000 — уже около 1-3%, и дальше погрешность уменьшается ещё больше.

    Для желающих самостоятельно оценить качество работы алгоритма — небольшой скриптик на perl-е:
    #!/usr/bin/perl -s
    
    $N = shift || 100000;
    for (1 .. $N) {
      push @a, (rand(1)*rand(1))**2;
    }
    @t = sort { $a <=> $b } @a;
    $MEDIAN = $t[$N/2-1];
    
    median_init();
    for (@a) {
      median_update($_);
    }
    
    $diff = ($MEDIAN-$median)/$MEDIAN*100;
    $avg = $sum/$count;
    
    print "Real:\t\t$MEDIAN
    Estimated:\t$median
    Difference:\t$diff \%
    
    Average: $avg
    ";
    
    sub median_init() {
      $count = $sum = $median = 0;
    }
    sub median_update($) {
      $value = $_[0];
      $sum += $value;
      $count ++;
      $delta = $sum / $count / $count;
      if($median <= $value) {
        $median += $delta;
        $s2 = '+';
      } else {
        $median -= $delta;
        $s2 = '-';
      }
      if($debug) {
        $s = ($value > $MEDIAN) ? '+' : '-';
        printf "%s %.5f\t%d\t%.7f\t  %s  %e\n", $s, $value, $count, $median, $s2, $delta;
      }
    }
    


    Статистика для разных распределений:

    rand(1) — медиана 0.5, среднее 0.5:
    размер выборки значения отклонений в процентах на 10-ти различных выборках
    100 -0.06568 -0.06159 -0.02009 -0.00372 -0.00177 -0.00012 0.00441 0.00503 0.01287 2.46360
    1000 -0.01196 -0.00462 -0.00415 -0.00240 0.00187 0.00948 0.01446 0.01787 0.02734 0.04150
    10000 -0.04025 -0.03712 -0.00983 -0.00512 0.00632 0.00768 0.01212 0.01267 0.05422 0.07043


    rand(1)^2 — медиана 0.25, среднее 1/3:
    размер 10 значений отклонений в процентах
    100 -0.80813 -0.38749 -0.38348 -0.32599 -0.31828 -0.13034 -0.12811 0.06157 0.28336 6.07339
    1000 -2.98383 -0.32917 -0.27234 -0.22337 -0.09080 0.03338 0.06851 0.30257 0.30286 0.31563
    10000 -0.97752 -0.51781 -0.44970 -0.39834 -0.17945 0.02517 0.06729 0.13353 0.31976 0.44967


    rand(1)*rand(1) — медиана 0.187.., среднее 0.25:
    размер 10 значений отклонений в процентах
    100 -3.84694 -0.08139 -0.06980 -0.04138 -0.02477 0.01441 0.02042 0.03998 0.16999 0.17690
    1000 -0.38108 -0.08865 -0.06287 -0.00716 0.01326 0.02373 0.05510 0.06785 0.09153 0.14421
    10000 -0.44392 -0.26293 -0.10081 -0.05271 0.02870 0.03629 0.09319 0.14522 0.14807 0.20980


    -log(rand(1)) — медиана 0.69..., среднее 1
    размер 10 значений отклонений в процентах
    100 -15.40835 -0.07608 -0.06993 -0.05958 -0.03355 -0.02186 -0.00486 0.01619 0.10614 0.11928
    1000 -3.05399 -0.06109 -0.05757 -0.04758 -0.01636 -0.01522 0.01376 0.03155 0.05036 0.06091
    10000 -0.24191 -0.10548 -0.09042 -0.04264 -0.03507 -0.01219 -0.01158 0.00481 0.04048 0.07969


    -log(rand(1)^4) — медиана 2.76.., среднее 4
    размер 10 значений отклонений в процентах
    100 -4.20326 -0.06970 -0.05433 -0.04191 -0.01110 -0.00353 0.02933 0.05837 0.06381 0.08264
    1000 -0.01198 -0.01090 -0.00719 0.01424 0.01454 0.01905 0.03147 0.03237 0.05578 0.89456
    10000 -0.45039 -0.06332 -0.02954 -0.02589 -0.01358 -0.01100 0.00229 0.01811 0.02675 0.05882
    Ads
    AdBlock has stolen the banner, but banners are not teeth — they will be back

    More

    Comments 22

      0
      Если нужна гарантированная точность 0.5%-2%, то я бы просто завел массив из 200 счетчиков и инкременировал бы нужный (координата считается одной операцией деления, которая у вас тоже есть). При выходе за границы — пересчет массива в 2 раза, в пределе очень редкое событие… 2 таких расширения в худшем случае дадут сужение в 4 раза. Дальше допилил бы алгоритм уже по ситуации.
      Или вы точность медианы не по крайним значениям считаете, а по сигме?
        +1
        Точность я считаю в процентах от действительного значения медианы (а не как отклонение вычисленного значения от середины выборки)

        С массивом не могли бы подробнее расписать? Мне кажется, Вы не учитываете некоторые моменты, которые возникают от того, что заранее неизвестно в каком диапазоне будут сосредоточены данные.
          +3
          Можно воспользоваться такой штукой github.com/HdrHistogram/HdrHistogram
          По факту это такой же массив, но шкала нелинейна. Можно хранить значения от 0 до 3,600,000,000 с точностью до трёх знаков.
          И считать не только медиану, но и другие квантили:)

          Есть порты на другие ЯП!;-)
            0
            >> в процентах от действительного значения медианы…
            Т.е. я правильно понимаю что если значения сосредоточены скажем в диапазоне от 100.0 до 101.0 то даже совсем криво посчитанная медиана (по крайнему значению) заведомо даст ошибку в не более чем 1% (например, если метод ошибочно дает 101.0 а реальная медиана это 100.01). И наоборот если реальная медиана топчется возле нуля — то даже небольшая абсолютная ошибка в относительном значении может болтаться от минус до плюс бесконечности.
            Но я просто не знаю реальной задачи, что там важно а что нет))

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

            Вообще, мне кажется точность медианы нужно оценивать по тому на сколько полученное разделение отличается от 50/50. Например требовать не хуже 49/51. Если так, то мне на вскидку кажется это невозможно выдержать не сохраняя и анализирую все прошлые значения ((
              +1
              На разных распределениях точность работы очень разная, Вы можете сами увидеть это из данных в конце статьи. Для положительных данных с длинным хвостом (схожих с распределением Ципфа) точность замечательная, но всегда можно придумать контрпример, где она будет ужасная. Другое дело, что встретить в реальности такой контрпример будет затруднительно.
            +1
            Нужно не только расширять интервал, но ещё и так же динамически увеличивать точность, добавлять счётчики — например, если диапазон от 0 до 1000, но большинство значений лежит между 0 и 1, то между 0 и 1 вам нужно соответственно больше счётчиков, иначе получите медиану где-то в районе ~5 вместо ~0,5.
              0
              Согласен. Но если точность считается от границ — то значение 5 вместо 0.5 приемлемо, т.к. это по прежнему — 0.5%. Выше написал дополнительно по этому поводу.
                0
                Точность медианы считать от границ диапазона?
                Не думаю, что в этом есть какой-то смысл.
            +2
            Просто оставлю тут вот это — ru.wikipedia.org/wiki/%D0%A1%D0%BA%D0%BE%D0%BB%D1%8C%D0%B7%D1%8F%D1%89%D0%B0%D1%8F_%D1%81%D1%80%D0%B5%D0%B4%D0%BD%D1%8F%D1%8F

            З. Ы. Добавьте, пожалуйста, таблички со сравнением вычисленной медианы по полной выборке, и вычисленной в итоге вашим методом на разных выборках.
              +1
              Добавил таблички, можете сравнить.
                0
                Я, может быть, чего-то не понял, но вроде бы автор говорил о медиане по всему набору данных, а в указанной статье речь в основном идёт об арифметическом среднем, причём скользящем, т.е. только по n последним значениям. Скользящие функции тривиально реализуются за O(n) времени и O(1) памяти, а вот по всему объёму — это уже интересней.
                  +1
                  На самом деле предложенный мною алгоритм тоже является в какой-то мере скользящим — результат больше зависит от последних значений, чем от первых.
                    0
                    Но первые ведь тоже учитываются, верно? Просто для меня это немного разный класс алгоритмов. Статистика по окну ничем не отличается от статистики по небольшому набору данных. А вот статистику по всему большому набору получить гораздо сложнее, просто в силу ограничений памяти и процессорного времени.
                      +1
                      Учитываются, верно. Просто потихоньку теряют свой вес.
                0
                Да, интересная задача. Если в один проход решать, то можно ткнуть в несколько случайных элементов и выбрать из них медиану. Если элементов будет порядка O(1/eps^2 \log(1/delta)), то вероятность, что выбранный элемент будет стоять от медианы на расстоянии не более eps * n будет не меньше 1 — delta.

                А еще вроде есть хитрый алгоритм Мунро-Патерсона, который умеет эту задачу решать точно с O(1) памяти за O(log n) проходов.
                  +2
                  Вопрос по теме и неплохой ответ на StackOverflow.
                    0
                    Да, ничто не ново под Луной ;) Я не особо усердно искал решение этой проблемы, ограничившись только рускоязычными материалами.
                    +1
                    Есть же Greenwald and Khanna GK01 и многочисленные его вариации, находятся по ключевым словам fast approximate quantiles.
                      0
                      По поводу эфективности… а оценка имеет наименьшую дисперсию в классе несмещённых оценок?

                      Такое ощущение, что вовсе нет.
                        +1
                        Уверен, что не имеет. Есть и более точные методы, но они требуют больше вычислений и ещё больше времени на имплементацию. Для того, чтобы глянуть, чему приблизительно равняется медиана, достаточно и такого простого метода.
                        0
                        А почему в одном случае из десяти так резко отличается отклонение? Ещё я не вижу сильной корреляции на разных размерах выборки, вы пробовали задавать размер в районе 10^6 — 10^8?
                          0
                          Случайное отклонение. На больших выборках (типа 10^6) их вероятность настолько ничтожна, что такие выбросы не встречаются.

                        Only users with full accounts can post comments. Log in, please.