Spitzensignalerfassung in Echtzeit-Zeitreihendaten

Question

---

Aktualisierung: Der bisher beste Algorithmus bisher ist dieser.

---

Diese Frage erkundet robuste Algorithmen zur Erkennung plötzlicher Peaks in echtzeitfähigen Zeitreihendaten.

Betrachten Sie das folgende Beispiel für Daten:

_{Beispiel dieser Daten ist im Matlab-Format (aber es handelt sich bei dieser Frage nicht um die Sprache, sondern um den Algorithmus):}

p = [1 1 1.1 1 0.9 1 1 1.1 1 0.9 1 1.1 1 1 0.9 1 1 1.1 1 1 1 1 1.1 0.9 1 1.1 1 1 0.9, ...
     1 1.1 1 1 1.1 1 0.8 0.9 1 1.2 0.9 1 1 1.1 1.2 1 1.5 1 3 2 5 3 2 1 1 1 0.9 1 1, ... 
     3 2.6 4 3 3.2 2 1 1 0.8 4 4 2 2.5 1 1 1];

Es ist deutlich zu erkennen, dass es drei große Peaks und einige kleine Peaks gibt. Diese Datensatz ist ein spezifisches Beispiel aus der Klasse von Zeitreihendatensätzen, um die es in der Frage geht. Diese Datenklasse hat zwei allgemeine Eigenschaften:

1. Es gibt grundlegendes Rauschen mit einem allgemeinen Mittelwert
2. Es gibt große 'Peaks' oder 'höhere Datenpunkte', die signifikant vom Rauschen abweichen.

Angenommen wird auch folgendes:

• Die Breite der Peaks kann nicht im Voraus bestimmt werden
• Die Höhe der Peaks weicht signifikant von den anderen Werten ab
• Der Algorithmus aktualisiert sich in Echtzeit (also mit jedem neuen Datenpunkt)

In einem solchen Szenario muss ein Grenzwert erstellt werden, der Signale auslöst. Allerdings kann der Grenzwert nicht statisch sein und muss in Echtzeit mithilfe eines Algorithmus bestimmt werden.

---

Meine Frage: Welcher Algorithmus ist gut geeignet, um solche Schwellenwerte in Echtzeit zu berechnen? Gibt es spezifische Algorithmen für solche Situationen? Was sind die bekanntesten Algorithmen?

---

_{Robuste Algorithmen oder nützliche Einsichten werden sehr geschätzt. (kann in jeder Sprache antworten: es geht um den Algorithmus)}

Answer 1

Hier ist eine C++-Implementierung des geglätteten Z-Score-Algorithmus aus dieser Antwort

std::vector smoothedZScore(std::vector input)
{   
    //Verzögerung 5 für die Glättungsfunktionen
    int lag = 5;
    //3,5 Standardabweichungen für das Signal
    float threshold = 3.5;
    //zwischen 0 und 1, wobei 1 normaler Einfluss ist, 0,5 ist die Hälfte
    float influence = .5;

    if (input.size() <= lag + 2)
    {
        std::vector emptyVec;
        return emptyVec;
    }

    //Variablen initialisieren
    std::vector signals(input.size(), 0.0);
    std::vector filteredY(input.size(), 0.0);
    std::vector avgFilter(input.size(), 0.0);
    std::vector stdFilter(input.size(), 0.0);
    std::vector subVecStart(input.begin(), input.begin() + lag);
    avgFilter[lag] = mean(subVecStart);
    stdFilter[lag] = stdDev(subVecStart);

    for (size_t i = lag + 1; i < input.size(); i++)
    {
        if (std::abs(input[i] - avgFilter[i - 1]) > threshold * stdFilter[i - 1])
        {
            if (input[i] > avgFilter[i - 1])
            {
                signals[i] = 1; //# Positives Signal
            }
            else
            {
                signals[i] = -1; //# Negatives Signal
            }
            //Machen Sie den Einfluss geringer
            filteredY[i] = influence* input[i] + (1 - influence) * filteredY[i - 1];
        }
        else
        {
            signals[i] = 0; //# Kein Signal
            filteredY[i] = input[i];
        }
        //Passen Sie die Filter an
        std::vector subVec(filteredY.begin() + i - lag, filteredY.begin() + i);
        avgFilter[i] = mean(subVec);
        stdFilter[i] = stdDev(subVec);
    }
    return signals;
}

Answer 2

Hier ist eine Implementierung des Smoothed z-score Algorithmus (oben) in Golang. Es geht von einem Slice von []int16 (PCM 16bit Samples) aus. Du kannst hier ein Gist finden.

/*
Einstellungen (die unten sind Beispiele: wähle das Beste für deine Daten)
setze Verzögerung auf 5;          # Verzögerung 5 für die Glättungsfunktionen
setze Schwellenwert auf 3,5;      # 3,5 Standardabweichungen für das Signal
setze Einfluss auf 0,5;           # zwischen 0 und 1, wobei 1 normaler Einfluss ist, 0,5 ist die Hälfte
*/

// ZScore auf 16bit WAV-Samples
func ZScore(samples []int16, lag int, threshold float64, influence float64) (signals []int16) {
    //lag := 20
    //threshold := 3.5
    //influence := 0.5

    signals = make([]int16, len(samples))
    filteredY := make([]int16, len(samples))
    for i, sample := range samples[0:lag] {
        filteredY[i] = sample
    }
    avgFilter := make([]int16, len(samples))
    stdFilter := make([]int16, len(samples))

    avgFilter[lag] = Average(samples[0:lag])
    stdFilter[lag] = Std(samples[0:lag])

    for i := lag + 1; i < len(samples); i++ {

        f := float64(samples[i])

        if float64(Abs(samples[i]-avgFilter[i-1])) > threshold*float64(stdFilter[i-1]) {
            if samples[i] > avgFilter[i-1] {
                signals[i] = 1
            } else {
                signals[i] = -1
            }
            filteredY[i] = int16(influence*f + (1-influence)*float64(filteredY[i-1]))
            avgFilter[i] = Average(filteredY[(i - lag):i])
            stdFilter[i] = Std(filteredY[(i - lag):i])
        } else {
            signals[i] = 0
            filteredY[i] = samples[i]
            avgFilter[i] = Average(filteredY[(i - lag):i])
            stdFilter[i] = Std(filteredY[(i - lag):i])
        }
    }

    return
}

// Durchschnittswert einer Menge von Werten
func Average(chunk []int16) (avg int16) {
    var sum int64
    for _, sample := range chunk {
        if sample < 0 {
            sample *= -1
        }
        sum += int64(sample)
    }
    return int16(sum / int64(len(chunk)))
}

Answer 3

C++ Implementierung

#include 
#include 
#include 
#include 
#include 
#include 
#include 

using namespace std;

typedef long double ld;
typedef unsigned int uint;
typedef std::vector::iterator vec_iter_ld;

/**
 * Überschreiben des ostream-Operators für die hübsche Formatierung von Vektoren.
 */
template
std::ostream &operator<<(std::ostream &os, std::vector vec) {
    os << "[";
    if (vec.size() != 0) {
        std::copy(vec.begin(), vec.end() - 1, std::ostream_iterator(os, " "));
        os << vec.back();
    }
    os << "]";
    return os;
}

/**
 * Diese Klasse berechnet den Mittelwert und die Standardabweichung eines Teilvektors.
 * Dies ist im Grunde die statistische Berechnung eines Teilvektors mit einer Fenstergröße "lag".
 */
class VektorStatistik {
public:
    /**
     * Konstruktor für die Klasse VektorStatistik.
     *
     * @param start - Dies ist die Iterator-Position des Starts des Fensters,
     * @param end   - Dies ist die Iterator-Position des Endes des Fensters,
     */
    VektorStatistik(vec_iter_ld start, vec_iter_ld end) {
        this->start = start;
        this->end = end;
        this->berechnen();
    }

    /**
     * Diese Methode berechnet den Mittelwert und die Standardabweichung mit der STL-Funktion.
     * Dies ist die Zwei-Pass-Implementierung der Berechnung von Mittelwert und Varianz.
     */
    void berechnen() {
        ld summe = std::accumulate(start, end, 0.0);
        uint schnittgröße = std::distance(start, end);
        ld mittelwert = summe / schnittgröße;
        std::vector diff(schnittgröße);
        std::transform(start, end, diff.begin(), [mittelwert](ld x) { return x - mittelwert; });
        ld sq_summe = std::inner_product(diff.begin(), diff.end(), diff.begin(), 0.0);
        ld standardabweichung = std::sqrt(sq_summe / schnittgröße);

        this->m1 = mittelwert;
        this->m2 = standardabweichung;
    }

    ld mittelwert() {
        return m1;
    }

    ld standardabweichung() {
        return m2;
    }

private:
    vec_iter_ld start;
    vec_iter_ld end;
    ld m1;
    ld m2;
};

/**
 * Dies ist die Implementierung des Smoothed Z-Score Algorithmus.
 * Dies ist die direkte Übersetzung von https://stackoverflow.com/a/22640362/1461896.
 *
 * @param input - Eingangssignal
 * @param lag - die Verzögerung des gleitenden Fensters
 * @param schwelle - der Z-Score, bei dem der Algorithmus ein Signal sendet
 * @param einfluss - der Einfluss (zwischen 0 und 1) neuer Signale auf den Mittelwert und die Standardabweichung
 * @return eine Map, die das gefilterte Signal und den entsprechenden Mittelwert und die Standardabweichung enthält.
 */
unordered_map> z_score_thresholding(vector input, int lag, ld schwelle, ld einfluss) {
    unordered_map> ausgabe;

    uint n = (uint) input.size();
    vector signale(input.size());
    vector gefilterter_input(input.begin(), input.end());
    vector gefilterter_mittelwert(input.size());
    vector gefilterte_stdabw(input.size());

    VektorStatistik lag_teilvektor_statistik(input.begin(), input.begin() + lag);
    gefilterter_mittelwert[lag - 1] = lag_teilvektor_statistik.mittelwert();
    gefilterte_stdabw[lag - 1] = lag_teilvektor_statistik.standardabweichung();

    for (int i = lag; i < n; i++) {
        if (abs(input[i] - gefilterter_mittelwert[i - 1]) > schwelle * gefilterte_stdabw[i - 1]) {
            signale[i] = (input[i] > gefilterter_mittelwert[i - 1]) ? 1.0 : -1.0;
            gefilterter_input[i] = einfluss * input[i] + (1 - einfluss) * gefilterter_input[i - 1];
        } else {
            signale[i] = 0.0;
            gefilterter_input[i] = input[i];
        }
        VektorStatistik lag_teilvektor_statistik(gefilterter_input.begin() + (i - lag), gefilterter_input.begin() + i);
        gefilterter_mittelwert[i] = lag_teilvektor_statistik.mittelwert();
        gefilterte_stdabw[i] = lag_teilvektor_statistik.standardabweichung();
    }

    ausgabe["signale"] = signale;
    ausgabe["gefilterter_mittelwert"] = gefilterter_mittelwert;
    ausgabe["gefilterte_stdabw"] = gefilterte_stdabw;

    return ausgabe;
};

int main() {
    vector input = {1.0, 1.0, 1.1, 1.0, 0.9, 1.0, 1.0, 1.1, 1.0, 0.9, 1.0, 1.1, 1.0, 1.0, 0.9, 1.0, 1.0, 1.1, 1.0,
                        1.0, 1.0, 1.0, 1.1, 0.9, 1.0, 1.1, 1.0, 1.0, 0.9, 1.0, 1.1, 1.0, 1.0, 1.1, 1.0, 0.8, 0.9, 1.0,
                        1.2, 0.9, 1.0, 1.0, 1.1, 1.2, 1.0, 1.5, 1.0, 3.0, 2.0, 5.0, 3.0, 2.0, 1.0, 1.0, 1.0, 0.9, 1.0,
                        1.0, 3.0, 2.6, 4.0, 3.0, 3.2, 2.0, 1.0, 1.0, 0.8, 4.0, 4.0, 2.0, 2.5, 1.0, 1.0, 1.0};

    int lag = 30;
    ld schwelle = 5.0;
    ld einfluss = 0.0;
    unordered_map> ausgabe = z_score_thresholding(input, lag, schwelle, einfluss);
    cout << ausgabe["signale"] << endl;
}

Answer 4

Ich dachte, ich würde meine Julia-Implementierung des Algorithmus für andere bereitstellen. Der Gist kann hier gefunden werden

using Statistics
using Plots
function SmoothedZscoreAlgo(y, lag, threshold, influence)
    # Julia-Implementierung von http://stackoverflow.com/a/22640362/6029703
    n = length(y)
    signals = zeros(n) # initiale Signalergebnisse
    filteredY = copy(y) # initiale gefilterte Serie
    avgFilter = zeros(n) # initiales Durchschnittsfilter
    stdFilter = zeros(n) # initiales Standardabweichungsfilter
    avgFilter[lag - 1] = mean(y[1:lag]) # initiale erste Wert
    stdFilter[lag - 1] = std(y[1:lag]) # initiale erste Wert

    for i in range(lag, stop=n-1)
        if abs(y[i] - avgFilter[i-1]) > threshold*stdFilter[i-1]
            if y[i] > avgFilter[i-1]
                signals[i] += 1 # positivess Signal
            else
                signals[i] += -1 # negatives Signal
            end
            # Mache den Einfluss kleiner
            filteredY[i] = influence*y[i] + (1-influence)*filteredY[i-1]
        else
            signals[i] = 0
            filteredY[i] = y[i]
        end
        avgFilter[i] = mean(filteredY[i-lag+1:i])
        stdFilter[i] = std(filteredY[i-lag+1:i])
    end
    return (signals = signals, avgFilter = avgFilter, stdFilter = stdFilter)
end

# Daten
y = [1,1,1.1,1,0.9,1,1,1.1,1,0.9,1,1.1,1,1,0.9,1,1,1.1,1,1,1,1,1.1,0.9,1,1.1,1,1,0.9,
       1,1.1,1,1,1.1,1,0.8,0.9,1,1.2,0.9,1,1,1.1,1.2,1,1.5,1,3,2,5,3,2,1,1,1,0.9,1,1,3,
       2.6,4,3,3.2,2,1,1,0.8,4,4,2,2.5,1,1,1]

# Einstellungen: lag = 30, threshold = 5, influence = 0
lag = 30
threshold = 5
influence = 0

results = SmoothedZscoreAlgo(y, lag, threshold, influence)
upper_bound = results[:avgFilter] + threshold * results[:stdFilter]
lower_bound = results[:avgFilter] - threshold * results[:stdFilter]
x = 1:length(y)

yplot = plot(x,y,color="blue", label="Y",legend=:topleft)
yplot = plot!(x,upper_bound, color="green", label="Obere Grenze",legend=:topleft)
yplot = plot!(x,results[:avgFilter], color="cyan", label="Durchschnittsfilter",legend=:topleft)
yplot = plot!(x,lower_bound, color="green", label="Untere Grenze",legend=:topleft)
signalplot = plot(x,results[:signals],color="red",label="Signale",legend=:topleft)
plot(yplot,signalplot,layout=(2,1),legend=:topleft)

Answer 5

Dieses Problem ähnelt einem, das ich in einem Hybrid-/Embedded-Systems-Kurs hatte, allerdings ging es dort um die Erkennung von Fehlern bei einem verrauschten Sensorsignal. Wir haben einen Kalman-Filter verwendet, um den versteckten Zustand des Systems zu schätzen/vorherzusagen, und dann statistische Analysen durchgeführt, um die Wahrscheinlichkeit für das Auftreten eines Fehlers zu bestimmen. Wir haben mit linearen Systemen gearbeitet, allerdings gibt es auch nichtlineare Varianten. Mir ist in Erinnerung geblieben, dass der Ansatz erstaunlich anpassungsfähig war, aber ein Modell der Systemdynamik erforderte.