Ermittlung von Mittelwert und Standardabweichung in Echtzeit

Was wäre der ideale Weg, um den Mittelwert und die Standardabweichung eines Signals für eine Echtzeitanwendung zu ermitteln? Ich möchte in der Lage sein, einen Controller auszulösen, wenn ein Signal für eine bestimmte Zeitspanne mehr als 3 Standardabweichungen vom Mittelwert aufweist.

Ich gehe davon aus, dass ein dedizierter DSP dies recht einfach tun würde, aber gibt es eine "Verknüpfung", die möglicherweise nicht so komplizierte Dinge erfordert?

Die Antwort von Jason R weist einen Fehler auf, der in Knuths "Art of Computer Programming" vol. 2. Das Problem tritt auf, wenn Sie eine Standardabweichung haben, die einen kleinen Bruchteil des Mittelwerts darstellt: Die Berechnung von E (x ^ 2) - (E (x) ^ 2) ist sehr empfindlich gegenüber Gleitkomma-Rundungsfehlern.

Sie können dies sogar selbst in einem Python-Skript versuchen:

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

Ich erhalte -128.0 als Antwort, was eindeutig nicht rechnerisch gültig ist, da die Mathematik vorhersagt, dass das Ergebnis nicht negativ sein sollte.

Knuth zitiert einen Ansatz (ich erinnere mich nicht an den Namen des Erfinders) zur Berechnung des laufenden Mittelwerts und der Standardabweichung, der ungefähr so ​​aussieht:

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

und dann ist der Wert von nach jedem Schritt mder Mittelwert, und die Standardabweichung kann als sqrt(S/n)oder sqrt(S/n-1)abhängig davon berechnet werden, welche Ihre bevorzugte Definition der Standardabweichung ist.

Die Gleichung, die ich oben schreibe, unterscheidet sich geringfügig von der in Knuth, ist jedoch rechnerisch äquivalent.

Wenn ich noch ein paar Minuten Zeit habe, codiere ich die obige Formel in Python und zeige, dass Sie eine nicht negative Antwort erhalten (die hoffentlich nahe am richtigen Wert liegt).

update: hier ist es.

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)

Ergebnis:

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}

Sie werden feststellen, dass es immer noch einige Rundungsfehler gibt, aber es ist nicht schlecht, während naive_statsnur Kotzen.

Bearbeiten: Gerade bemerkt, dass Belisarius Wikipedia zitiert, in dem der Knuth-Algorithmus erwähnt wird.

Was wäre der ideale Weg, um den Mittelwert und die Standardabweichung eines Signals für eine Echtzeitanwendung zu ermitteln? Ich möchte in der Lage sein, einen Controller auszulösen, wenn ein Signal für eine bestimmte Zeitspanne mehr als 3 Standardabweichungen vom Mittelwert aufweist.

Der richtige Ansatz in solchen Situationen besteht normalerweise darin, einen exponentiell gewichteten laufenden Durchschnitt und eine Standardabweichung zu berechnen. Im exponentiell gewichteten Durchschnitt werden die Schätzungen des Mittelwerts und der Varianz in Richtung der letzten Stichprobe verschoben, wodurch Sie Schätzungen des Mittelwerts und der Varianz über die letzten Sekunden erhalten. $\tau$ ist wahrscheinlich das, was Sie möchten, und nicht der übliche arithmetische Durchschnitt über alle Stichproben jemals gesehen.

Im Frequenzbereich ist ein "exponentiell gewichteter laufender Durchschnitt" einfach ein echter Pol. Es ist einfach im Zeitbereich zu implementieren.

Zeitbereichsimplementierung

Sei meanund meansqsei die aktuelle Schätzung des Mittelwerts und des Mittelwerts des Quadrats des Signals. Aktualisieren Sie diese Schätzungen in jedem Zyklus mit der neuen Stichprobe 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);

Hier ist eine Konstante, die die effektive Länge des laufenden Durchschnitts bestimmt. So wählen Sie eine unten in „Analyse“ beschrieben wird.$0 < a < 1$$a$

Was oben als Imperativprogramm ausgedrückt wurde, kann auch als Signalflussdiagramm dargestellt werden:

Analyse

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

Wenn Sie die IIR-Filter in eigenen Blöcken zusammenfassen, sieht das Diagramm nun folgendermaßen aus:

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

$a$

Verweise

• Die Simulink-Diagrammquelle kann von https://gist.github.com/1942771 heruntergeladen werden

Eine Methode, die ich zuvor in einer eingebetteten Verarbeitungsanwendung verwendet habe, besteht darin, Akkumulatoren der Summe und der Quadratsumme des interessierenden Signals zu verwalten:

$i$$i$

oder Sie können verwenden:

je nachdem, welche Methode zur Schätzung der Standardabweichung Sie bevorzugen . Diese Gleichungen basieren auf der Definition der Varianz :

Ich habe diese in der Vergangenheit erfolgreich verwendet (obwohl ich mich nur mit der Varianzschätzung befasst habe, nicht mit der Standardabweichung), obwohl Sie vorsichtig mit den numerischen Typen sein müssen, die Sie zum Halten der Akkumulatoren verwenden, wenn Sie eine Summierung durchführen möchten eine lange Zeitspanne; du willst nicht überlaufen.

Bearbeiten: Zusätzlich zu dem obigen Kommentar zum Überlauf sollte beachtet werden, dass dies kein numerisch robuster Algorithmus ist, wenn er in Gleitkomma-Arithmetik implementiert ist, was möglicherweise große Fehler in der geschätzten Statistik verursacht. Schauen Sie sich die Antwort von Jason S an, um einen besseren Ansatz in diesem Fall zu finden.

Ähnlich wie bei der obigen bevorzugten Antwort (Jason S.) und auch abgeleitet von der von Knut übernommenen Formel (Vol.2, S. 232) kann man auch eine Formel ableiten, um einen Wert zu ersetzen, dh einen Wert in einem Schritt zu entfernen und hinzuzufügen . Nach meinen Tests liefert "Ersetzen" eine bessere Präzision als die zweistufige Version "Entfernen / Hinzufügen".

Der folgende Code ist in Java, meanund saktualisiert werden ( „global“ Membervariablen), gleich wie mund soben in Jasons Post. Der Wert countbezieht sich auf die Fenstergröße 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;
}

Die Antwort von Jason und Nibot unterscheidet sich in einem wichtigen Aspekt: ​​Die Methode von Jason berechnet den Standardwert und den Mittelwert für das gesamte Signal (da y = 0), während die von Nibot eine "laufende" Berechnung ist, dh, sie wiegt neuere Abtastwerte, die stärker sind als Abtastwerte aus der ferne Vergangenheit.

Da die Anwendung std dev und mean als Funktion der Zeit zu erfordern scheint, ist die Methode von Nibot wahrscheinlich die geeignetere (für diese spezielle Anwendung). Der wirklich schwierige Teil wird jedoch darin bestehen, den Teil mit der richtigen Zeitgewichtung zu finden. Das Beispiel von Nibot verwendet einen einfachen einpoligen Filter.

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

Die Auswahl des Tiefpassfilters richtet sich nach dem, was Sie über Ihr Signal wissen und nach der Zeitauflösung, die Sie für Ihre Schätzung benötigen. Niedrigere Grenzfrequenzen und eine höhere Ordnung führen zu einer besseren Genauigkeit, aber einer langsameren Reaktionszeit.

Um die Sache noch weiter zu komplizieren, wird ein Filter im linearen Bereich und ein anderer im quadratischen Bereich angewendet. Durch Quadrieren wird der spektrale Inhalt des Signals erheblich geändert, sodass Sie möglicherweise einen anderen Filter im Quadratbereich verwenden möchten.

Hier ist ein Beispiel, wie man Mittelwert, Effektivwert und Standardabweichung als Funktion der Zeit abschätzt.

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