Określanie średniej i odchylenia standardowego w czasie rzeczywistym

Jaki byłby idealny sposób na znalezienie średniej i odchylenia standardowego sygnału dla aplikacji w czasie rzeczywistym. Chciałbym być w stanie wyzwolić kontroler, gdy sygnał był większy niż 3 odchylenie standardowe od średniej przez pewien czas.

Zakładam, że dedykowany DSP zrobiłby to dość łatwo, ale czy jest jakiś „skrót”, który może nie wymagać czegoś tak skomplikowanego?

W odpowiedzi Jasona R jest błąd, który jest omawiany w „Art of Computer Programming” Knutha, t. 2. Problem pojawia się, gdy masz odchylenie standardowe, które stanowi niewielki ułamek średniej: obliczenie E (x ^ 2) - (E (x) ^ 2) cierpi z powodu dużej wrażliwości na błędy zaokrąglania zmiennoprzecinkowego.

Możesz nawet spróbować tego samodzielnie w skrypcie 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

Otrzymuję -128,0 jako odpowiedź, co oczywiście nie jest poprawne obliczeniowo, ponieważ matematyka przewiduje, że wynik powinien być nieujemny.

Knuth przytacza podejście (nie pamiętam nazwiska wynalazcy) do obliczania średniej bieżącej i odchylenia standardowego, które wygląda mniej więcej tak:

 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);

a następnie po każdym kroku wartość mjest średnią, a odchylenie standardowe można obliczyć jako sqrt(S/n)lub w sqrt(S/n-1)zależności od ulubionej definicji odchylenia standardowego.

Równanie, które piszę powyżej, jest nieco inne niż w Knuth, ale jest obliczeniowo równoważne.

Kiedy mam jeszcze kilka minut, koduję powyższą formułę w Pythonie i pokażę, że otrzymasz nieujemną odpowiedź (która, mam nadzieję, jest zbliżona do poprawnej wartości).

aktualizacja: oto jest.

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)

wynik:

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}

Zauważysz, że nadal występuje błąd zaokrąglania, ale nie jest zły, podczas gdy naive_statspo prostu wymiotuje.

edycja: Właśnie zauważyłem komentarz Belizariusza cytujący Wikipedię, która wspomina o algorytmie Knuth.

Jaki byłby idealny sposób na znalezienie średniej i odchylenia standardowego sygnału dla aplikacji w czasie rzeczywistym. Chciałbym być w stanie wyzwolić kontroler, gdy sygnał był większy niż 3 odchylenie standardowe od średniej przez pewien czas.

$\tau$

W dziedzinie częstotliwości „wykładniczo ważona średnia bieżąca” jest po prostu prawdziwym biegunem. Jest prosty do wdrożenia w dziedzinie czasu.

Implementacja w dziedzinie czasu

Niech meani meansqbyć obecne szacunki średniej i średniej kwadratu sygnału. W każdym cyklu aktualizuj te szacunki o nową próbkę 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$

To, co zostało wyrażone powyżej jako program imperatywny, może być również przedstawione jako diagram przepływu sygnału:

Analiza

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

Scalając filtry IIR we własne bloki, schemat wygląda teraz tak:

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

$a$

Referencje

• Źródło diagramu Simulink można pobrać ze strony https://gist.github.com/1942771

Metodą, której użyłem wcześniej w aplikacji do przetwarzania wbudowanego, jest utrzymanie akumulatorów sumy i sumy kwadratów sygnału zainteresowania:

$i$$i$

lub możesz użyć:

w zależności od preferowanej metody szacowania odchylenia standardowego . Te równania są oparte na definicji wariancji :

Używałem ich z powodzeniem w przeszłości (chociaż byłem zainteresowany jedynie oszacowaniem wariancji, a nie odchyleniem standardowym), chociaż musisz uważać na typy numeryczne, których używasz do przechowywania akumulatorów, jeśli zamierzasz sumować długi okres czasu; nie chcesz przepełnienia.

Edycja: Oprócz powyższego komentarza dotyczącego przepełnienia należy zauważyć, że nie jest to algorytm odporny numerycznie, gdy jest implementowany w arytmetyki zmiennoprzecinkowej, potencjalnie powodując duże błędy w szacowanych statystykach. Spójrz na odpowiedź Jasona S, aby uzyskać lepsze podejście w tej sprawie.

Podobnie do powyższej preferowanej odpowiedzi (Jason S.), a także wywodzącej się ze wzoru zaczerpniętego z Knuta (Vol. 2, str. 232), można również wyprowadzić formułę, aby zastąpić wartość, tj. Usunąć i dodać wartość w jednym kroku . Według moich testów, replace zapewnia lepszą precyzję niż dwuetapowa wersja „usuń / dodaj”.

Poniższy kod jest w Javie meani jest saktualizowany („globalne” zmienne składowe), tak samo jak mi spowyżej w poście Jasona. Wartość countodnosi się do rozmiaru okna 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;
}

Odpowiedź Jasona i Nibota różni się jednym ważnym aspektem: metoda Jasona oblicza std dev i średnią dla całego sygnału (od y = 0), podczas gdy Nibota jest obliczeniem „bieżącym”, tj. Waży nowsze próbki silniejsze niż próbki z odległa przeszłość.

Ponieważ wydaje się, że aplikacja wymaga std dev i średniej w funkcji czasu, metoda Nibota jest prawdopodobnie bardziej odpowiednia (dla tej konkretnej aplikacji). Jednak prawdziwą trudną częścią będzie prawidłowe wyważenie czasu. Przykład Nibota wykorzystuje prosty filtr jednobiegunowy.

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

Wybór filtra dolnoprzepustowego może zależeć od tego, co wiesz o swoim sygnale i rozdzielczości czasu potrzebnej do oszacowania. Niższe częstotliwości odcięcia i wyższy porządek będą skutkować lepszą dokładnością, ale wolniejszym czasem reakcji.

Aby jeszcze bardziej skomplikować sprawę, jeden filtr stosuje się w dziedzinie liniowej, a drugi w dziedzinie do kwadratu. Kwadratowanie znacząco zmienia zawartość spektralną sygnału, więc możesz chcieć użyć innego filtra w domenie kwadratu.

Oto przykład, jak oszacować średnią, rms i std dev w funkcji czasu.

%% 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);