Как стать автором
Обновить

Комментарии 12

Спасибо за подробный разбор алгоритма! Но вообще-то почти в точности такая же технология используется уже несколько десятилетий, и называется "обработка в скользящем окне". Например, мы у себя (при работе с временными рядами) используем подобный подход с 1992 года. При добавлении нового элемента в выборку все суммы инкрементируются, при удалении элемента - декрементируются. При этом как поступающий в выборку, так и удаляемый из нее элемент в нашем случае может быть так называемым пропуском, который кодируется Nan или другим особым значением. Поэтому перед инкрементом/декрементом сумм выполняется дополнительная проверка на Nan.

А вот деревья мы не используем, так как абсолютные моменты бывают нужны

очень редко,

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

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

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

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

из паспорта сигнала

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

Это позволяет избежать вычитания двух больших близких чисел, и очень существенно повышает точность.

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

Полные тексты статей

которые упомянуты выше, можно найти вот тут.

а сама программа

в виде exe-шников доступна для скачивания вот здесь,

А ее исходники и файлы проектов - вот тут.

Спасибо за ваш интерес!

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

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

Кстати, по поводу вычислительной стабильности, работа 2016-го года, на которую мы ссылаемся в тексте (https://link.springer.com/article/10.1007/s00180-015-0637-z) о быстром обновлении центральных моментов (не абсолютных) почти целиком посвящена проблемам накопления ошибки вычислений с плавающей точкой. Этот вопрос там хорошо и подробно рассмотрен, рекомендуем.

А откуда берутся значения с плавающей точкой? Если речь о данных из "реального мира", то обычно они приходят с АЦП, т.е. по факту целочисленны. В этом случае можно все накопление/обновление сумм выполнять в целочисленной арифметике и никаких проблем с точностью не будет.

А откуда берутся значения с плавающей точкой? 

Это от задачи зависит.

Например, у нас в геофизическом мониторинге основная задача - это изучение природных явлений и закономерностей, а не

управление ими

То есть, речь чаще всего идет о поиске каких-то эффектов или проверке гипотез по ретроспективным данным. Обычно для этого анализируются сигналы, накопленные за много лет. А они копятся в физических единицах, естественно - так как разряды АЦП после каждой перенастройки аппаратуры разные, а ряд должен быть однородным. Не говоря уже о том, что для хранения коэффициентов перевода (если в базе лежит выходной сигнал АЦП) потребуется отдельная нехилая БД...

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

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

И в таких задачах без человека не обойтись. А везде, где участвует человек, гораздо удобнее работать в физических единицах - градусах там или нанострейнах, а не в условных "отсчетах АЦП"

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

Соответственно, данные для всяких расчетов берутся не с АЦП, а из БД, где они лежат уже в виде значений с плавающей точкой.

Но понятно, что случаи бывают разные. Рутинные задачи реального времени очень часто допускают работу в integer. И тогда этим надо пользоваться, разумеется

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

Ну даже если данные реально с плавающей точкой, но диапазон их значений известен, то возможно проще их квантовать и считать статистику по квантованным целочисленным значениям

Это понятная схема, но она не всегда применима на практике, так как плотность real-шкалы неравномерна, в отличие от integer. И очень часто это существенно. Так что проблема вовсе не сводится к ошибке квантования (в обычном ее понимании).

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

еще на предварительной стадии (до записи в базу).

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

под каким углом мы его изучаем

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

Иногда для идентификации выбросов достаточно тонкие алгоритмы приходится применять.

И еще один ключевой момент.

Если Вы продумали, что в результате хитрых трюков первичной обработки у нас в итоге все равно сохраняется примерно однозначное соответствие между исходными int-значениями с выхода АЦП и "физическими" real-значениями, то

это совершенно не так

При геофизическом мониторинге длина рядов может достигать десятков лет. За это время аппаратура и установка могут меняться неоднократно. А уж замена/регулировка усилителя (АЦП) - это и вовсе рутинная процедура. Например, когда наклон в результате дрейфа уплывает за границу диапазона, мы меняем смещение (электронно или механически), возвращая сигнал обратно к центру шкалы. После этого все формулы пересчета в физ. единицы меняются, и тот же самый отсчет АЦП будет соответствовать совершенно другой физической величине. Другая типичная ситуация - это коррекция коэффициента усиления. Мы же не знаем заранее, как себя поведет сигнал с уникального датчика в уникальных условиях. Поэтому время от времени может возникнуть необходимость расширить динамический диапазон (если сигналу его не хватает) или, наоборот, сузить (если сигнал упорно занимает лишь небольшую часть полной шкалы). Это тоже приводит к пересчету всех формул.

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

мы этого не делаем

И причина тут даже не одна, а все две:

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

2) Ну и во-вторых, мы свои алгоритмы модифицируем иногда. И сделать это гораздо проще, когда алгоритм существует в одном варианте, чем когда существует зоопарк версий...

Но, разумеется, если бы мы делали узкоспециализированный блок для обработки сигналов строго определенного типа с заранее известными свойствами, и вопрос бы уперся в скорость обработки, то такую оптимизацию следовало бы рассмотреть в числе первых

В общем, спасибо Вам за желание подкинуть конструктивный совет - кому-то он здесь наверняка пригодится (ведь случайные люди такую статью не полезут читать ;-).

Но вот именно нам - увы...

Таким образом, у нас получилось обновить k-й абсолютный центральный момент выборки размера n за время O(k\log n), хотя для этого нам пришлось содержать сбалансированное дерево с n вершинами и с аккумулятором размера O(k) в каждой вершине.

Кажется, есть метод попроще, без всяких кастомных деревьев. Достаточно даже встроенного в язык heap/set. Да даже если их нет, написать кучу гораздо проще и константа у нее будет поменьше.

Идея в том, чтобы поддерживать 2 дерева/множества - одно для чисел меньше-равных среднего, второе - для чисел больше. Все операции разбиваются на 2 части: сначала добавляем новые числа и пересчитываем среднее, потом изменяем деревья, чтобы числа были на своих местах.

Добавление числа совсем просто - сравниваем со средним и кладем или туда или сюда. Пересчет среднего - обновление одного аккумулятора. Потом "балансировка" - вынимаем из левого дерева максимум, пока он больше среднего и кладем в правое дерево. И наоборот: из правого дерева достаем минимум. Да, тут должно стать понятно, что у вас два heap с разными знаками - левый ищет максимум, правый - минимум.

O(k) аккумуляторов тут надо всего 2 - для левого и правого множества целиком. Пересчитываем их при перемещении элементов между множествами.

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

Только может ттак получиться, что много элементов каждый раз будут менять heap. Например, если данные вот такие: 1, 0, 0, -1.1, 1.1, 0, 0, -1.1, 1.1, .... то среднее всё время скачет от "чуть больше нуля" до "чуть меньше нуля" и все нули нужно перекладывать между хипами.

Да, действительно.

Формула бинома в статье приведена. Ее небольшое упрощение для n=2 - нет, потому что оно не обобщается на абсолютные моменты.

извините - проглядел

Ее небольшое упрощение для n=2 

Там всегда последние пару слагаемых упрощаются в формуле (если вы не заметили)

...+ (-1)^(n-1)*С(n-1,n)*Ex*(Ex)^(n-1)+ (-1)^n*C(n,n)*(Ex)^n = ... + ...

оно не обобщается на абсолютные моменты.

))))) Вы правы - только в половине случаев: n=2*k => E(x-Ex)^n = E|x-Ex|^n && E(x^n)=E(|x|^n)

Зарегистрируйтесь на Хабре, чтобы оставить комментарий