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

31

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

Я предполагаю, что выделенный DSP сделает это довольно легко, но есть ли какой-нибудь «ярлык», который не требует чего-то столь сложного?

statistics real-time measurement jonsca
источник

Вы знаете что-нибудь о сигнале? Это стационарно?

@Tim Давайте скажем, что это стационарно. Для моего собственного любопытства, каковы будут последствия нестационарного сигнала?

Джонса

3

Если он стационарный, вы можете просто вычислить среднее значение и стандартное отклонение. Все было бы сложнее, если бы среднее и стандартное отклонение менялось со временем.

5

Очень похожие: en.wikipedia.org/wiki/…

доктор Белизариус

36

В ответе Джейсона Р. есть недостаток, который обсуждается в книге Кнута «Искусство компьютерного программирования», том. 2. Проблема возникает, если у вас есть стандартное отклонение, которое составляет небольшую долю от среднего значения: вычисление E (x ^ 2) - (E (x) ^ 2) страдает от серьезной чувствительности к ошибкам округления с плавающей запятой.

Вы даже можете попробовать это самостоятельно в скрипте Python:

ofs = 1e9
A = [ofs+x for x in [1,-1,2,3,0,4.02,5]] 
A2 = [x*x for x in A]
(sum(A2)/len(A))-(sum(A)/len(A))**2

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

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

 initialize:
    m = 0;
    S = 0;
    n = 0;

 for each incoming sample x:
    prev_mean = m;
    n = n + 1;
    m = m + (x-m)/n;
    S = S + (x-m)*(x-prev_mean);

и затем после каждого шага значение mявляется средним, и стандартное отклонение может быть рассчитано как sqrt(S/n)или в sqrt(S/n-1)зависимости от того, какое ваше любимое определение стандартного отклонения.

Уравнение, которое я пишу выше, немного отличается от уравнения Кнута, но оно эквивалентно в вычислительном отношении.

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

Обновление: вот оно.

test1.py:

import math

def stats(x):
  n = 0
  S = 0.0
  m = 0.0
  for x_i in x:
    n = n + 1
    m_prev = m
    m = m + (x_i - m) / n
    S = S + (x_i - m) * (x_i - m_prev)
  return {'mean': m, 'variance': S/n}

def naive_stats(x):
  S1 = sum(x)
  n = len(x)
  S2 = sum([x_i**2 for x_i in x])
  return {'mean': S1/n, 'variance': (S2/n - (S1/n)**2) }

x1 = [1,-1,2,3,0,4.02,5] 
x2 = [x+1e9 for x in x1]

print "naive_stats:"
print naive_stats(x1)
print naive_stats(x2)

print "stats:"
print stats(x1)
print stats(x2)

результат:

naive_stats:
{'variance': 4.0114775510204073, 'mean': 2.0028571428571427}
{'variance': -128.0, 'mean': 1000000002.0028572}
stats:
{'variance': 4.0114775510204073, 'mean': 2.0028571428571431}
{'variance': 4.0114775868357446, 'mean': 1000000002.0028571}

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

редактировать: только что заметил комментарий Велисария со ссылкой на Википедию, в которой упоминается алгоритм Кнута.

Джейсон С
источник

1

+1 за подробный ответ с примером кода. Этот подход превосходит тот, который указан в моем ответе, когда требуется реализация с плавающей запятой.

Джейсон Р

1

Можно также проверить это для реализации C ++: johndcook.com/standard_deviation.html

Руи Маркиз

1

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

Джейсон С

«Кнут цитирует подход (я не помню имя изобретателя) для расчета среднего значения» - кстати, это метод Уэлфорда .

Джейсон С

Если кто-то может помочь, я отправил вопрос, связанный с этим: dsp.stackexchange.com/questions/31812/…

Джонатан

13

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

$\tau$

В частотной области «экспоненциально взвешенное скользящее среднее» является просто реальным полюсом. Это просто реализовать во временной области.

Реализация во временной области

Позвольте meanи meansqбыть текущие оценки среднего и среднего квадрата сигнала. На каждом цикле обновляйте эти оценки новым образцом x:

% update the estimate of the mean and the mean square:
mean = (1-a)*mean + a*x
meansq = (1-a)*meansq + a*(x^2)

% calculate the estimate of the variance:
var = meansq - mean^2;

% and, if you want standard deviation:
std = sqrt(var);

$0 < a < 1$ $a$

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

введите описание изображения здесь

Анализ

$y_i = a x_i + (1-a) y_{i-1}$ $x_i$ $i$ $y_i$ $z$

H (z) = \frac{a}{1 - (1 - a) z^{- 1}}

$H(z) = \frac{a}{1-(1-a)z^{-1}}$

Сжимая БИХ-фильтры в их собственные блоки, диаграмма теперь выглядит так:

введите описание изображения здесь

$z = e^{s T}$ $T$ $f_s = 1/T$ $1-(1-a)e^{-sT}=0$ $s = \frac{1}{T} \log (1-a)$

$a$

a = 1 - \exp {2 π \frac{T}{τ}}

$a = 1 - \exp \left\{2\pi\frac{T}{\tau}\right\}$

Ссылки

Исходный код Simulink можно скачать с https://gist.github.com/1942771.

nibot
источник

1

a

$a$

a

$a$ 0 > a > 1

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

Шнарф

Woops! Конечно я имел ввиду 0 < a < 1. Если в вашей системе есть выборка tmie, Tи вы хотите усреднить постоянную времени tau, выберите a = 1 - exp (2*pi*T/tau).

нибот

Я думаю, что здесь может быть ошибка. Однополюсные фильтры не имеют коэффициента усиления 0 дБ при постоянном токе, и поскольку вы применяете один фильтр в линейной области и один в квадрате, ошибка усиления для E <x> и E <x ^ 2> различна. Я уточню в своем ответе

Хильмар

Он имеет усиление 0 дБ при постоянном токе. Замените z=1(DC) на H(z) = a/(1-(1-a)/z)и вы получите 1.

nibot

5

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

A_{x, i} = \sum_{k = 0}^{i} x [k] = A_{x, i - 1} + x [i], A_{x, - 1} = 0

$A_{x,i} = \sum_{k=0}^{i}x[k] = A_{x,i-1} + x[i], A_{x,-1} = 0$

A_{x^{2}, i} = \sum_{k = 0}^{i} x^{2} [k] = A_{x^{2}, i - 1} + x^{2} [i], A_{x^{2}, - 1} = 0

$A_{x^2,i} = \sum_{k=0}^{i}x^2[k] = A_{x^2,i-1} + x^2[i], A_{x^2,-1} = 0$

$i$ $i$

\tilde{μ} = \frac{A_{x_{i}}}{i + 1}

$\tilde\mu = \frac{A_{x_i}}{i+1}$

\tilde{σ} = \sqrt{\frac{A_{x_{i}^{2}}}{i + 1} - {\tilde{μ}}^{2}}

$\tilde\sigma = \sqrt{\frac{A_{x^2_i}}{i+1} - \tilde\mu^2}$

или вы можете использовать:

\tilde{σ} = \sqrt{\frac{A_{x_{i}^{2}}}{i} - {\tilde{μ}}^{2}}

$\tilde\sigma = \sqrt{\frac{A_{x^2_i}}{i} - \tilde\mu^2}$

в зависимости от того, какой метод оценки стандартного отклонения вы предпочитаете . Эти уравнения основаны на определении дисперсии :

σ^{2} = E (X^{2}) - (E (X))^{2}

$\sigma^2 = \operatorname{E}(X^2) - (\operatorname{E}(X))^2$

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

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

Джейсон Р
источник

1

A_{x, i} = x [i] + A_{x, i - 1}, A_{x, 0} = x [0]

$A_{x,i}=x[i]+A_{x,i-1},\ A_{x,0}=x[0]$

i

$i$

x

$x$

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

Джейсон Р

2

-1, когда у меня достаточно репутации, чтобы сделать это: у этого есть числовые проблемы. Смотри Кнут том. 2

Джейсон С

σ

$\sigma$

μ^{2}

$\mu^2$

σ^{2} = E (X^{2}) - (E (X))^{2}

$\sigma^2 = E(X^2) - (E(X))^2$

2

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

Джейсон Р

3

Подобно предпочтительному ответу выше (Jason S.), а также полученному из формулы, взятой из Кнута (Vol.2, p 232), можно также получить формулу для замены значения, то есть удалить и добавить значение за один шаг. , Согласно моим тестам, замена обеспечивает лучшую точность, чем двухэтапная версия удаления / добавления.

Код ниже написан на Java meanи sобновляется («глобальные» переменные-члены) так же, как mи sвыше в посте Джейсона. Значение countотносится к размеру окна n.

/**
 * Replaces the value {@code x} currently present in this sample with the
 * new value {@code y}. In a sliding window, {@code x} is the value that
 * drops out and {@code y} is the new value entering the window. The sample
 * count remains constant with this operation.
 * 
 * @param x
 *            the value to remove
 * @param y
 *            the value to add
 */
public void replace(double x, double y) {
    final double deltaYX = y - x;
    final double deltaX = x - mean;
    final double deltaY = y - mean;
    mean = mean + deltaYX / count;
    final double deltaYp = y - mean;
    final double countMinus1 = count - 1;
    s = s - count / countMinus1 * (deltaX * deltaX - deltaY * deltaYp) - deltaYX * deltaYp / countMinus1;
}

марко
источник

3

Ответ Джейсона и Нибота отличается в одном важном аспекте: метод Джейсона вычисляет стандартное отклонение и среднее значение для всего сигнала (так как y = 0), в то время как ответ Нибота является «бегущим» вычислением, то есть он взвешивает более свежие выборки сильнее, чем выборки из далекое прошлое.

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

$E[x]$ $x[n]$ $E[x^2]$ $x[n]^2$

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

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

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

%% example
fs = 44100; n = fs; % 44.1 kHz sample rate, 1 second
% signal: white noise plus a low frequency drift at 5 Hz)
x = randn(n,1) + sin(2*pi*(0:n-1)'*5/fs);
% mean estimation filter: since we are looking for effects in the 5 Hz range we use maybe a
% 25 Hz filter, 2nd order so it's not too sluggish
[b,a] = butter(2,25*2/fs);
xmeanEst = filter(b,a,x);
% now we estimate x^2, since most frequency double we use twice the bandwidth
[b,a] = butter(2,50*2/fs);
x2Est = filter(b,a,x.^2);
% std deviation estimate
xstd = sqrt(x2Est)-xmeanEst;
% and plot it
h = plot([x, xmeanEst sqrt(x2Est) xstd]);
grid on;
legend('x','E<x>','sqrt(E<x^2>)','Std dev');
set(h(2:4),'Linewidth',2);

Hilmar
источник

1

Фильтр в моем ответе соответствует y1 = filter(a,[1 (1-a)],x);.

нибот

1

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

Джейсон Р

nibot, прости, что ты прав, а я был неправ. Я исправлю это прямо сейчас

Хильмар

1

+1 за предложение другой фильтрации для х и х ^ 2

нибот

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

Ответы: