Онлайн алгоритм для среднего абсолютного отклонения и большого набора данных

У меня есть небольшая проблема, которая заставляет меня волноваться. Я должен написать процедуру для онлайн-процесса приобретения многомерного временного ряда. На каждом временном интервале (например, 1 секунда) я получаю новую выборку, которая в основном представляет собой вектор с плавающей запятой размера N. Операция, которую мне нужно сделать, немного сложнее:

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

2. Я рассчитываю средний процент вектора таким же образом, но с использованием прошлых значений.

3. Для каждого прошедшего значения я вычисляю абсолютное отклонение вектора процентов, относящегося к этой выборке, с вектором глобального среднего процента, вычисленным на шаге 2. Таким образом, абсолютное отклонение всегда равно 0 (когда вектор равен среднему вектор) и 2 (когда он полностью отличается).

4. Используя среднее значение отклонений для всех предыдущих выборок, я вычисляю среднее абсолютное отклонение, которое снова является числом от 0 до 2.

5. Я использую среднее абсолютное отклонение, чтобы определить, совместим ли новый образец с другими образцами (сравнивая его абсолютное отклонение со средним абсолютным отклонением всего набора, вычисленным на шаге 4).

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

Спасибо за вашу помощь.

Если вы можете принять некоторую неточность, эта проблема может быть легко решена путем подсчета биннинга . То есть выберите какое-то большое число (скажем, M = 1000 ), затем инициализируйте несколько целочисленных бинов B i , j для i = 1 … M и j = 1 … N , где N - размер вектора, равный нулю. Затем, когда вы увидите k- е наблюдение процентного вектора, увеличьте B i , j, если - 1 ) /$M$$M = 1000$$B_{i,j}$$i = 1\ldots M$$j = 1\ldots N$$N$$k$$B_{i,j}$ й элемент этого вектора находится между$j$$(i-1)/M$ и , перебирая N элементов вектора. (Я предполагаю, что ваши входные векторы неотрицательны, поэтому, когда вы вычисляете свои «проценты», векторы находятся в диапазоне [ 0 , 1 ] .)$i/M$$N$$[0,1]$

В любой момент времени вы можете оценить средний вектор из бункеров и среднее абсолютное отклонение. После наблюдения таких векторов, то J - й элемент средней оценивается ˉ Х J = $K$$j$иJй элемент среднего абсолютного отклонения оценивается

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

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

ВНИМАНИЕ: Это не онлайн алгоритм. Это требует O(n)памяти. Кроме того, он имеет худшую производительность O(n), например, для наборов данных [1, -2, 4, -8, 16, -32, ...](т. Е. Такой же, как и полный пересчет). [1]

Тем не менее, поскольку он по-прежнему хорошо работает во многих случаях, его стоит опубликовать здесь. Например, чтобы рассчитать абсолютное отклонение 10000 случайных чисел от -100 до 100 при поступлении каждого элемента, мой алгоритм занимает менее одной секунды, в то время как полный пересчет занимает более 17 секунд (на моем компьютере будет зависеть от машины и согласно входным данным). Однако вам необходимо сохранить весь вектор в памяти, что может быть ограничением для некоторых применений. Схема алгоритма следующая:

1. Вместо одного вектора для хранения прошлых измерений используйте три отсортированные приоритетные очереди (что-то вроде кучи мин / макс). Эти три списка делят входные данные на три: элементы больше среднего, элементы меньше среднего и элементы равные среднему.
2. (Почти) каждый раз, когда вы добавляете элемент, среднее значение меняется, поэтому нам нужно переделить. Важнейшей вещью является отсортированный характер разделов, что означает, что вместо сканирования каждого элемента в списке для повторного использования нам нужно только читать те элементы, которые мы перемещаем. Хотя в худшем случае это все равно потребуетO(n) операций перемещения, для многих случаев это не так.
3. Используя некоторую умную бухгалтерию, мы можем убедиться, что отклонение правильно рассчитано во все времена, при перераспределении и при добавлении новых элементов.

Пример кода на Python приведен ниже. Обратите внимание, что он позволяет только элементы, которые будут добавлены в список, но не удалены. Это можно легко добавить, но в то время, когда я писал это, мне это не нужно. Вместо того, чтобы реализовывать очереди приоритетов самостоятельно, я использовал отсортированный список из превосходного пакета blist Даниэля Штутцбаха , который использует B + Tree для внутреннего использования.

Рассмотрим этот код под лицензией MIT . Он не был значительно оптимизирован или отшлифован, но работал для меня в прошлом. Новые версии будут доступны здесь . Дайте мне знать, если у вас есть какие-либо вопросы, или найдите какие-либо ошибки.

from blist import sortedlist
import operator

class deviance_list:
    def __init__(self):
        self.mean =  0.0
        self._old_mean = 0.0
        self._sum =  0L
        self._n =  0  #n items
        # items greater than the mean
        self._toplist =  sortedlist()
        # items less than the mean
        self._bottomlist = sortedlist(key = operator.neg)
        # Since all items in the "eq list" have the same value (self.mean) we don't need
        # to maintain an eq list, only a count
        self._eqlistlen = 0

        self._top_deviance =  0
        self._bottom_deviance =  0

    @property
    def absolute_deviance(self):
        return self._top_deviance + self._bottom_deviance

    def append(self,  n):
        # Update summary stats
        self._sum += n
        self._n +=  1
        self._old_mean =  self.mean
        self.mean =  self._sum /  float(self._n)

        # Move existing things around
        going_up = self.mean > self._old_mean
        self._rebalance(going_up)

        # Add new item to appropriate list
        if n >  self.mean:
            self._toplist.add(n)
            self._top_deviance +=  n -  self.mean
        elif n == self.mean: 
            self._eqlistlen += 1
        else:
            self._bottomlist.add(n)
            self._bottom_deviance += self.mean -  n


    def _move_eqs(self,  going_up):
        if going_up:
            self._bottomlist.update([self._old_mean] *  self._eqlistlen)
            self._bottom_deviance += (self.mean - self._old_mean) * self._eqlistlen
            self._eqlistlen = 0
        else:
            self._toplist.update([self._old_mean] *  self._eqlistlen)
            self._top_deviance += (self._old_mean - self.mean) * self._eqlistlen
            self._eqlistlen = 0


    def _rebalance(self, going_up):
        move_count,  eq_move_count = 0, 0
        if going_up:
            # increase the bottom deviance of the items already in the bottomlist
            if self.mean !=  self._old_mean:
                self._bottom_deviance += len(self._bottomlist) *  (self.mean -  self._old_mean)
                self._move_eqs(going_up)


            # transfer items from top to bottom (or eq) list, and change the deviances
            for n in iter(self._toplist):
                if n < self.mean:
                    self._top_deviance -= n -  self._old_mean
                    self._bottom_deviance += (self.mean -  n)
                    # we increment movecount and move them after the list
                    # has finished iterating so we don't modify the list during iteration
                    move_count +=  1
                elif n == self.mean:
                    self._top_deviance -= n -  self._old_mean
                    self._eqlistlen += 1
                    eq_move_count +=  1
                else:
                    break
            for _ in xrange(0,  move_count):
                self._bottomlist.add(self._toplist.pop(0))
            for _ in xrange(0,  eq_move_count):
                self._toplist.pop(0)

            # decrease the top deviance of the items remain in the toplist
            self._top_deviance -= len(self._toplist) *  (self.mean -  self._old_mean)
        else:
            if self.mean !=  self._old_mean:
                self._top_deviance += len(self._toplist) *  (self._old_mean -  self.mean)
                self._move_eqs(going_up)
            for n in iter(self._bottomlist): 
                if n > self.mean:
                    self._bottom_deviance -= self._old_mean -  n
                    self._top_deviance += n -  self.mean
                    move_count += 1
                elif n == self.mean:
                    self._bottom_deviance -= self._old_mean -  n
                    self._eqlistlen += 1
                    eq_move_count +=  1
                else:
                    break
            for _ in xrange(0,  move_count):
                    self._toplist.add(self._bottomlist.pop(0))
            for _ in xrange(0,  eq_move_count):
                self._bottomlist.pop(0)

            # decrease the bottom deviance of the items remain in the bottomlist
            self._bottom_deviance -= len(self._bottomlist) *  (self._old_mean -  self.mean)


if __name__ ==  "__main__":
    import random
    dv =  deviance_list()
    # Test against some random data,  and calculate result manually (nb. slowly) to ensure correctness
    rands = [random.randint(-100,  100) for _ in range(0,  1000)]
    ns = []
    for n in rands: 
        dv.append(n)
        ns.append(n)
        print("added:%4d,  mean:%3.2f,  oldmean:%3.2f,  mean ad:%3.2f" %
              (n, dv.mean,  dv._old_mean,  dv.absolute_deviance / dv.mean))
        assert sum(ns) == dv._sum,  "Sums not equal!"
        assert len(ns) == dv._n,  "Counts not equal!"
        m = sum(ns) / float(len(ns))
        assert m == dv.mean,  "Means not equal!"
        real_abs_dev = sum([abs(m - x) for x in ns])
        # Due to floating point imprecision, we check if the difference between the
        # two ways of calculating the asb. dev. is small rather than checking equality
        assert abs(real_abs_dev - dv.absolute_deviance) < 0.01, (
            "Absolute deviances not equal. Real:%.2f,  calc:%.2f" %  (real_abs_dev,  dv.absolute_deviance))

[1] Если симптомы не проходят, обратитесь к врачу.

Существует также параметрический подход. Игнорируя векторную природу ваших данных и рассматривая только маргиналы, достаточно решить проблему: найти онлайновый алгоритм для вычисления среднего абсолютного отклонения скаляра$X$, Если (и это большое «если» здесь) вы думали, что$X$Следуя некоторому распределению вероятности с неизвестными параметрами, вы можете оценить параметры, используя некоторый онлайн-алгоритм, а затем вычислить среднее абсолютное отклонение на основе этого параметризованного распределения. Например, если вы думали, что$X$ было (приблизительно) нормально распределено, вы можете оценить его стандартное отклонение, как $s$и среднее абсолютное отклонение будет оцениваться как $s \sqrt{2 / \pi}$(см. Половинное нормальное распределение ).

MAD (x) - это всего лишь два одновременных вычисления медианы, каждое из которых может быть выполнено в режиме онлайн с помощью binmedian алгоритма .

Вы можете найти соответствующую статью, а также код C и FORTRAN онлайн здесь .

(это всего лишь использование хитрого трюка поверх хитрого трюка Шаббишефа, чтобы сэкономить память).

Приложение:

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

Нижеследующее обеспечивает неточное приближение, хотя неточность будет зависеть от распределения входных данных. Это онлайн-алгоритм, но он только приближает абсолютное отклонение. Он основан на хорошо известном алгоритме для расчета дисперсии онлайн, описанном Уэлфордом в 1960-х годах. Его алгоритм, переведенный в R, выглядит так:

M2 <- 0
mean <- 0
n <- 0

var.online <- function(x){
    n <<- n + 1
    diff <- x - mean
    mean <<- mean + diff / n
    M2 <<- M2 + diff * (x - mean)
    variance <- M2 / (n - 1)
    return(variance)
}

Он работает очень похоже на встроенную функцию дисперсии R:

set.seed(2099)
n.testitems <- 1000
n.tests <- 100
differences <- rep(NA, n.tests)
for (i in 1:n.tests){
        # Reset counters
        M2 <- 0
        mean <- 0
        n <- 0

        xs <- rnorm(n.testitems)
        for (j in 1:n.testitems){
                v <- var.online(xs[j])
        }

        differences[i] <- abs(v - var(xs))

}
summary(differences)
     Min.   1st Qu.    Median      Mean   3rd Qu.      Max. 
0.000e+00 2.220e-16 4.996e-16 6.595e-16 9.992e-16 1.887e-15 

Модификация алгоритма для вычисления абсолютного отклонения просто включает дополнительный sqrtвызов. Тем не менее, sqrtвводит неточности, которые отражаются в результате:

absolute.deviance.online <- function(x){
    n <<- n + 1
    diff <- x - mean
    mean <<- mean + diff / n
    a.dev <<- a.dev + sqrt(diff * (x - mean))
    return(a.dev)
}

Ошибки, рассчитанные, как указано выше, намного больше, чем для расчета дисперсии:

    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
0.005126 0.364600 0.808000 0.958800 1.360000 3.312000 

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