Określanie średniej i odchylenia standardowego w czasie rzeczywistym

31

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?

jonsca
źródło
Czy wiesz coś o tym sygnale? Czy to jest stacjonarne?
@ Tim Powiedzmy, że jest stacjonarny. Co do mojej własnej ciekawości, jakie byłyby konsekwencje niestacjonarnego sygnału?
jonsca
3
Jeśli jest stacjonarny, możesz po prostu obliczyć średnią bieżącą i odchylenie standardowe. Sprawa byłaby bardziej skomplikowana, gdyby średnia i odchylenie standardowe zmieniały się z czasem.
5
Bardzo spokrewniony: en.wikipedia.org/wiki/…
Dr. belisarius

Odpowiedzi:

36

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.

Jason S.
źródło
1
+1 za szczegółową odpowiedź z przykładowym kodem. Podejście to jest lepsze od wskazanego w mojej odpowiedzi, gdy potrzebna jest implementacja zmiennoprzecinkowa.
Jason R
1
Można to również sprawdzić pod kątem implementacji w C ++: johndcook.com/standard_deviation.html
Rui Marques
1
tak, to wszystko. Używa dokładnych równań, które stosuje Knuth. Możesz nieco zoptymalizować i uniknąć konieczności sprawdzania początkowej iteracji w porównaniu do kolejnych iteracji, jeśli używasz mojej metody.
Jason S
„Knuth cytuje podejście (nie pamiętam nazwiska wynalazcy) do obliczania średniej biegu” - tak przy okazji, to metoda Welforda .
Jason S
Zadałem
Jonathan
13

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.

τ

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<1a

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

wprowadź opis zdjęcia tutaj

Analiza

yi=axi+(1a)yi1xiiyiz

H(z)=a1(1a)z1

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

wprowadź opis zdjęcia tutaj

z=esTTfs=1/T1(1a)esT=0s=1Tlog(1a)

a

a=1exp{2πTτ}

Referencje

nibot
źródło
1
aa0 > a > 1
Jest to podobne do podejścia Jasona R. Ta metoda będzie mniej dokładna, ale trochę szybsza i mniej pamięci. Takie podejście kończy się zastosowaniem okna wykładniczego.
sznurek
Woops! Oczywiście, że miałem na myśli 0 < a < 1. Jeśli twój system ma próbkowanie Tmie Ti chciałbyś uśrednić stałą czasową tau, wybierz a = 1 - exp (2*pi*T/tau).
nibot
Myślę, że może tu być błąd. Filtry jednobiegunowe nie mają wzmocnienia 0 dB przy DC, a ponieważ stosuje się jeden filtr w dziedzinie liniowej i jeden w dziedzinie kwadratu, błąd wzmocnienia jest inny dla E <x> i E <x ^ 2>. Rozwinę więcej w mojej odpowiedzi
Hilmar
Ma wzmocnienie 0 dB przy DC. Zastępować z=1(DC) do H(z) = a/(1-(1-a)/z)i masz 1.
nibot
5

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:

Ax,i=k=0ix[k]=Ax,i1+x[i],Ax,1=0

Ax2,i=k=0ix2[k]=Ax2,i1+x2[i],Ax2,1=0

ii

μ~=Axii+1

σ~=Axi2i+1μ~2

lub możesz użyć:

σ~=Axi2iμ~2

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

σ2=E(X2)(E(X))2

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.

Jason R.
źródło
1
Ax,i=x[i]+Ax,i1, Ax,0=x[0]ix
Tak, tak jest lepiej. Próbowałem przepisać, aby bardziej przejrzysta była implementacja rekurencyjna.
Jason R
2
-1, gdy mam wystarczającą liczbę powtórzeń: ma to problemy numeryczne. Zobacz Knuth vol. 2
Jason S
σμ2σ2=E(X2)(E(X))2
2
@JasonS: Nie zgadzam się z tym, że technika jest z natury wadliwa, chociaż zgadzam się z twoim twierdzeniem, że nie jest to metoda niezawodna numerycznie, gdy jest stosowana w zmiennoprzecinkowym. Powinienem był powiedzieć, że z powodzeniem wykorzystałem to w aplikacji, która używa arytmetyki liczb całkowitych . Arytmetyka liczb całkowitych (lub stacjonarnych implementacji liczb ułamkowych) nie cierpi z powodu wskazanego przez ciebie problemu, który powoduje utratę precyzji. W tym kontekście jest to odpowiednia metoda, która wymaga mniej operacji na próbkę.
Jason R
3

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;
}
marco
źródło
3

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[x2]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);
Hilmar
źródło
1
Filtr w mojej odpowiedzi odpowiada y1 = filter(a,[1 (1-a)],x);.
nibot
1
Dobra uwaga na temat rozróżnienia między bieżącymi statystykami a statystykami z ogólnej próby. Moją implementację można zmodyfikować w celu obliczania statystyk bieżących poprzez kumulowanie się nad ruchomym oknem, co można również zrobić skutecznie (na każdym etapie odejmuj próbkę czasu, która właśnie wysunęła się z okna z każdego akumulatora).
Jason R
nibot, przepraszam, masz rację, a ja się myliłem. Zaraz to poprawię
Hilmar,
1
+1 za sugerowanie innego filtrowania dla xi x ^ 2
nibot