Алгоритм отклонения подвижного состава

Ответ 1

Я столкнулся с этой проблемой. Есть несколько отличных сообщений в вычислениях кумулятивной дисперсии, таких как John Cooke. Точно вычисляет сообщение об ошибке и сообщение из Digital explorations, код Python для вычисления выборок и дисперсий населения, ковариации и коэффициента корреляции. Просто не удалось найти никаких адаптированных к качению окна.

"Выполнение стандартных отклонений" от почтовых сообщений было критическим в том, чтобы заставить формулу качения оконтировать. Джим берет энергетическую сумму квадратов разностей значений по сравнению с методом Уолфорда с использованием суммы квадратов различий среднего. Формула выглядит следующим образом:

PSA today = PSA (вчера) + (((x сегодня * x сегодня) - x вчера))/n

• x = значение в ваших временных рядах
• n = количество проанализированных вами значений.

Но для преобразования формулы Power Sum Average в оконный сорт вам нужно подстроить формулу следующим образом:

PSA today = PSA вчера + (((x сегодня * x сегодня) - (x вчера * x вчера)/n

• x = значение в ваших временных рядах
• n = количество проанализированных вами значений.

Вам также понадобится формула Rolling Simple Moving Average:

SMA today = SMA вчера + ((x сегодня - x сегодня - n)/n

• x = значение в ваших временных рядах
• n = период, используемый для вашего окна перехода.

Оттуда вы можете вычислить отклонение опрокидывания:

Население Var сегодня = (сегодня PSA * n - n * SMA сегодня * SMA)/n

Или отклонение от опрокидывания:

Пример Var сегодня = (PSA сегодня * n - n * SMA сегодня * SMA сегодня)/(n - 1)

Я рассмотрел эту тему вместе с образцом кода Python в сообщении в блоге несколько лет назад, Running Variance.

Надеюсь, что это поможет.

Обратите внимание: я предоставил ссылки на все сообщения в блогах и математические формулы в латекс (изображения) для этого ответа. Но из-за моей низкой репутации (< 10); Я ограничен только двумя гиперссылками и абсолютно без изображений. сожалею об этом. Надеюсь, это не уберет контент.

Ответ 2

Я занимаюсь той же проблемой.

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

next_index = (index + 1) % window_size;    // oldest x value is at next_index, wrapping if necessary.

new_mean = mean + (x_new - xs[next_index])/window_size;

Я адаптировал алгоритм Welford и работает для всех значений, с которыми я тестировал.

varSum = var_sum + (x_new - mean) * (x_new - new_mean) - (xs[next_index] - mean) * (xs[next_index] - new_mean);

xs[next_index] = x_new;
index = next_index;

Чтобы получить текущую дисперсию, просто разделите varSum на размер окна: variance = varSum / window_size;

Ответ 3

Если вы предпочитаете использовать код над словами (в основном на основе записи DanS): http://calcandstuff.blogspot.se/2014/02/rolling-variance-calculation.html

public IEnumerable RollingSampleVariance(IEnumerable data, int sampleSize)
{
    double mean = 0;
    double accVar = 0;

    int n = 0;
    var queue = new Queue(sampleSize);

    foreach(var observation in data)
    {
        queue.Enqueue(observation);
        if (n < sampleSize)
        {
            // Calculating first variance
            n++;
            double delta = observation - mean;
            mean += delta / n;
            accVar += delta * (observation - mean);
        }
        else
        {
            // Adjusting variance
            double then = queue.Dequeue();
            double prevMean = mean;
            mean += (observation - then) / sampleSize;
            accVar += (observation - prevMean) * (observation - mean) - (then - prevMean) * (then - mean);
        }

        if (n == sampleSize)
            yield return accVar / (sampleSize - 1);
    }
}

Ответ 4

Здесь применяется подход с разделителем и покорением, который имеет O(log k) -time обновления, где k - количество выборок. Он должен быть относительно стабильным по тем же причинам, что парное суммирование и БПФ стабильны, но это немного сложнее, а константа невелика.

Предположим, что мы имеем последовательность A длины m со средним E(A) и дисперсией V(A), а последовательность B длины n со средним значением E(B) и дисперсией V(B). Пусть C - конкатенация A и B. Мы имеем

p = m / (m + n)
q = n / (m + n)
E(C) = p * E(A) + q * E(B)
V(C) = p * (V(A) + (E(A) + E(C)) * (E(A) - E(C))) + q * (V(B) + (E(B) + E(C)) * (E(B) - E(C)))

Теперь нарисуйте элементы в красно-черном дереве, где каждый node украшен средним значением и дисперсией поддерева, внедренного в этот node. Вставить справа; удалите слева. (Поскольку мы получаем доступ только к концам, дерево splay может быть O(1) амортизировано, но я предполагаю, что амортизация является проблемой для вашего приложения.) Если k известно во время компиляции, возможно, вы можете развернуть внутренний цикл FFTW.

Ответ 5

На самом деле алгоритм Weldords позволяет AFAICT легко адаптироваться для вычисления взвешенной дисперсии. И, установив весы на -1, вы должны иметь возможность эффективно отменить элементы. Я не проверял математику, разрешает ли она отрицательные веса, но при первом взгляде это должно быть!

Я сделал небольшой эксперимент, используя ELKI:

void testSlidingWindowVariance() {
MeanVariance mv = new MeanVariance(); // ELKI implementation of weighted Welford!
MeanVariance mc = new MeanVariance(); // Control.

Random r = new Random();
double[] data = new double[1000];
for (int i = 0; i < data.length; i++) {
  data[i] = r.nextDouble();
}

// Pre-roll:
for (int i = 0; i < 10; i++) {
  mv.put(data[i]);
}
// Compare to window approach
for (int i = 10; i < data.length; i++) {
  mv.put(data[i-10], -1.); // Remove
  mv.put(data[i]);
  mc.reset(); // Reset statistics
  for (int j = i - 9; j <= i; j++) {
    mc.put(data[j]);
  }
  assertEquals("Variance does not agree.", mv.getSampleVariance(),
    mc.getSampleVariance(), 1e-14);
}
}

Я получаю около 14 цифр точности по сравнению с точным двухпроходным алгоритмом; это примерно столько же, сколько можно ожидать от парных. Обратите внимание, что Welford действительно получает некоторые вычислительные затраты из-за дополнительных делений - он занимает примерно в два раза больше, чем точный двухпроходный алгоритм. Если размер вашего окна невелик, может быть гораздо более разумным фактически пересчитать среднее значение, а затем во втором пропускать дисперсию каждый раз.

Я добавил этот эксперимент как unit test в ELKI, вы можете увидеть полный источник здесь: http://elki.dbs.ifi.lmu.de/browser/elki/trunk/test/de/lmu/ifi/dbs/elki/math/TestSlidingVariance.java он также сравнивается с точной двухпроходной дисперсией.

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

Ответ 6

Я с нетерпением жду того, что я ошибаюсь, но я не думаю, что это можно сделать "быстро". Тем не менее, большая часть расчета отслеживает EV над окном, что можно сделать легко.

Я уйду с вопросом: вы уверены, что вам нужна оконная функция? Если вы не работаете с очень большими окнами, вероятно, лучше всего использовать известный заранее определенный алгоритм.

Ответ 7

Я предполагаю, что отслеживаю ваши 20 образцов, Sum (X ^ 2 от 1..20) и Sum (X от 1..20), а затем последовательно пересчитывая две суммы на каждой итерации, недостаточно эффективен? Можно повторно пересчитывать новую дисперсию без сложения, возведения в квадрат и т.д. Всех образцов каждый раз.

Как в:

Sum(X^2 from 2..21) = Sum(X^2 from 1..20) - X_1^2 + X_21^2
Sum(X from 2..21) = Sum(X from 1..20) - X_1 + X_21

Ответ 8

Здесь другое решение O(log k): найдите квадраты исходной последовательности, затем пары сумм, затем четверки и т.д. (вам нужно немного буфера, чтобы можно было найти все это эффективно.) Затем добавьте те значения, которые вам нужны, чтобы получить ответ. Например:

|||||||||||||||||||||||||  // Squares
| | | | | | | | | | | | |  // Sum of squares for pairs
|   |   |   |   |   |   |  // Pairs of pairs
|       |       |       |  // (etc.)
|               |
   ^------------------^    // Want these 20, which you can get with
        |       |          // one...
    |   |       |   |      // two, three...
                    | |    // four...
   ||                      // five stored values.

Теперь вы используете стандартную формулу E (x ^ 2) -E (x) ^ 2, и все готово. (Нет, если вам нужна хорошая стабильность для небольших наборов чисел, это было предполагая, что это было только скопление ошибок при прокачке, которые вызывали проблемы.)

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

Ответ 9

Только для 20 значений, тривиально, чтобы адаптировать метод, открытый здесь (я не сказал быстро, хотя).

Вы можете просто выбрать массив из 20 из этих классов RunningStat.

Первые 20 элементов потока несколько особенные, но как только это делается, это намного проще:

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

Очевидно, что этот подход не является масштабируемым...

Вы также можете заметить, что в числах, которые мы сохраняем, есть некоторая redudancy (если вы идете с полным классом RunningStat). Очевидным улучшением было бы сохранить 20 последних Mk и Sk напрямую.

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

Ответ 10

Я знаю, что этот вопрос старый, но в случае, если кто-то заинтересован в этом, следует код Python. Он вдохновлен johndcook в блоге, @Joachim, @DanS и комментарии @Jaime. Приведенный ниже код все еще дает небольшие неточности для небольших размеров окон данных. Наслаждайтесь.

from __future__ import division
import collections
import math


class RunningStats:
    def __init__(self, WIN_SIZE=20):
        self.n = 0
        self.mean = 0
        self.run_var = 0
        self.WIN_SIZE = WIN_SIZE

        self.windows = collections.deque(maxlen=WIN_SIZE)

    def clear(self):
        self.n = 0
        self.windows.clear()

    def push(self, x):

        self.windows.append(x)

        if self.n <= self.WIN_SIZE:
            # Calculating first variance
            self.n += 1
            delta = x - self.mean
            self.mean += delta / self.n
            self.run_var += delta * (x - self.mean)
        else:
            # Adjusting variance
            x_removed = self.windows.popleft()
            old_m = self.mean
            self.mean += (x - x_removed) / self.WIN_SIZE
            self.run_var += (x + x_removed - old_m - self.mean) * (x - x_removed)

    def get_mean(self):
        return self.mean if self.n else 0.0

    def get_var(self):
        return self.run_var / (self.WIN_SIZE - 1) if self.n > 1 else 0.0

    def get_std(self):
        return math.sqrt(self.get_var())

    def get_all(self):
        return list(self.windows)

    def __str__(self):
        return "Current window values: {}".format(list(self.windows))