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

Wir haben versucht, den geglätteten Z-Score-Algorithmus auf unseren Datensatz anzuwenden, was entweder zu einer Überempfindlichkeit oder Unterempfindlichkeit führt (abhängig davon, wie die Parameter eingestellt sind), mit wenig Mittelweg. In unserem Website-Verkehrssignal haben wir eine niedrige Frequenz-Baseline beobachtet, die den täglichen Zyklus darstellt und selbst mit den bestmöglichen Parametern (unten dargestellt) hat sie immer noch besonders am 4. Tag nachgelassen, weil die meisten Datenpunkte als Anomalie erkannt werden.

Basierend auf dem ursprünglichen Z-Score-Algorithmus haben wir eine Möglichkeit entwickelt, dieses Problem durch Rückwärtsfilterung zu lösen. Einzelheiten des modifizierten Algorithmus und dessen Anwendung auf die Zuordnung von Fernsehwerbetraffic sind auf unserem Team-Blog veröffentlicht.

Answer 2

Anhang 1 zur Originalantwort: Matlab und R Übersetzungen

Matlab-Code

function [signale,avgFilter,stdFilter] = ThresholdingAlgo(y,verzögerung,schwelle,einfluss)
% Initialisiere Signale
signale = zeros(length(y),1);
% Initialisiere gefilterte Serie
gefiltertY = y(1:verzögerung+1);
% Initialisiere Filter
avgFilter(verzögerung+1,1) = mean(y(1:verzögerung+1));
stdFilter(verzögerung+1,1) = std(y(1:verzögerung+1));
% Schleife über alle Datenpunkte y(verzögerung+2),...,y(t)
for i=verzögerung+2:length(y)
    % Wenn neuer Wert eine bestimmte Anzahl von Abweichungen entfernt ist
    if abs(y(i)-avgFilter(i-1)) > schwelle*stdFilter(i-1)
        if y(i) > avgFilter(i-1)
            % Positives Signal
            signale(i) = 1;
        else
            % Negatives Signal
            signale(i) = -1;
        end
        % Mache Einfluss geringer
        gefiltertY(i) = einfluss*y(i)+(1-einfluss)*gefiltertY(i-1);
    else
        % Kein Signal
        signale(i) = 0;
        gefiltertY(i) = y(i);
    end
    % Passe die Filter an
    avgFilter(i) = mean(gefiltertY(i-verzögerung:i));
    stdFilter(i) = std(gefiltertY(i-verzögerung:i));
end
% Fertig, nun Ergebnisse zurückgeben
end

Beispiel:

% 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
verzögerung = 30;
schwelle = 5;
einfluss = 0;

% Ergebnisse erhalten
[signale,avg,dev] = ThresholdingAlgo(y,verzögerung,schwelle,einfluss);

figure; subplot(2,1,1); hold on;
x = 1:length(y); ix = verzögerung+1:length(y);
area(x(ix),avg(ix)+schwelle*dev(ix),'FaceColor',[0.9 0.9 0.9],'EdgeColor','none');
area(x(ix),avg(ix)-schwelle*dev(ix),'FaceColor',[1 1 1],'EdgeColor','none');
plot(x(ix),avg(ix),'LineWidth',1,'Color','cyan','LineWidth',1.5);
plot(x(ix),avg(ix)+schwelle*dev(ix),'LineWidth',1,'Color','green','LineWidth',1.5);
plot(x(ix),avg(ix)-schwelle*dev(ix),'LineWidth',1,'Color','green','LineWidth',1.5);
plot(1:length(y),y,'b');
subplot(2,1,2);
stairs(signale,'r','LineWidth',1.5); ylim([-1.5 1.5);

---

R-Code

ThresholdingAlgo <- function(y,verzögerung,schwelle,einfluss) {
  signale <- rep(0,length(y))
  gefiltertY <- y[1:verzögerung]
  avgFilter <- NULL
  stdFilter <- NULL
  avgFilter[verzögerung] <- mean(y[1:verzögerung], na.rm=TRUE)
  stdFilter[verzögerung] <- sd(y[1:verzögerung], na.rm=TRUE)
  for (i in (verzögerung+1):length(y)){
    if (abs(y[i]-avgFilter[i-1]) > schwelle*stdFilter[i-1]) {
      if (y[i] > avgFilter[i-1]) {
        signale[i] <- 1;
      } else {
        signale[i] <- -1;
      }
      gefiltertY[i] <- einfluss*y[i]+(1-einfluss)*gefiltertY[i-1]
    } else {
      signale[i] <- 0
      gefiltertY[i] <- y[i]
    }
    avgFilter[i] <- mean(gefiltertY[(i-verzögerung):i], na.rm=TRUE)
    stdFilter[i] <- sd(gefiltertY[(i-verzögerung):i], na.rm=TRUE)
  }
  return(list("signale"=signale,"avgFilter"=avgFilter,"stdFilter"=stdFilter))
}

Beispiel:

# Daten
y <- c(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)

verzögerung <- 30
schwelle <- 5
einfluss <- 0

# Algorithmus mit verzögerung = 30, schwelle = 5, einfluss = 0 ausführen
ergebnis <- ThresholdingAlgo(y,verzögerung,schwelle,einfluss)

# Ergebnis plotten
par(mfrow = c(2,1),oma = c(2,2,0,0) + 0.1,mar = c(0,0,2,1) + 0.2)
plot(1:length(y),y,type="l",ylab="",xlab="") 
lines(1:length(y),ergebnis$avgFilter,type="l",col="cyan",lwd=2)
lines(1:length(y),ergebnis$avgFilter+schwelle*ergebnis$stdFilter,type="l",col="green",lwd=2)
lines(1:length(y),ergebnis$avgFilter-schwelle*ergebnis$stdFilter,type="l",col="green",lwd=2)
plot(ergebnis$signale,type="S",col="red",ylab="",xlab="",ylim=c(-1.5,1.5),lwd=2)

Dieser Code (beide Sprachen) liefert für die Daten der Originalfrage das folgende Ergebnis:

---

Anhang 2 zur Originalantwort: Matlab Demo-Code

(Klicken Sie, um Daten zu erstellen)

function [] = RobustThresholdingDemo()

%% SPEZIFIKATIONEN
verzögerung = 5;       % Verzögerung für die Glättung
schwelle   = 3.5;     % Anzahl der Standardabweichungen vom Mittelwert, um ein Signal zu erkennen
einfluss   = 0.3;     % Bei Signal: Wie viel Einfluss für neue Daten? (zwischen 0 und 1)
                       % 1 ist normaler Einfluss, 0,5 ist halb      
%% STARTEN SIE DAS DEMO
DemoScreen(30,verzögerung,schwelle,einfluss);

end

function [signale,avgFilter,stdFilter] = ThresholdingAlgo(y,verzögerung,schwelle,einfluss)
signale = zeros(length(y),1);
gefiltertY = y(1:verzögerung+1);
avgFilter(verzögerung+1,1) = mean(y(1:verzögerung+1));
stdFilter(verzögerung+1,1) = std(y(1:verzögerung+1));
for i=verzögerung+2:length(y)
    if abs(y(i)-avgFilter(i-1)) > schwelle*stdFilter(i-1)
        if y(i) > avgFilter(i-1)
            signale(i) = 1;
        else
            signale(i) = -1;
        end
        gefiltertY(i) = einfluss*y(i)+(1-einfluss)*gefiltertY(i-1);
    else
        signale(i) = 0;
        gefiltertY(i) = y(i);
    end
    avgFilter(i) = mean(gefiltertY(i-verzögerung:i));
    stdFilter(i) = std(gefiltertY(i-verzögerung:i));
end
end

% Demo-Bildschirmfunktion
function [] = DemoScreen(n,verzögerung,schwelle,einfluss)
figure('Position',[200 100,1000,500]);
subplot(2,1,1);
title(sprintf(['Zeichnen von Datenpunkten (%.0f max)      [Einstellungen: Verzögerung = %.0f, '...
    'Schwelle = %.2f, Einfluss = %.2f]'],n,verzögerung,schwelle,einfluss));
ylim([0 5]); xlim([0 50]);
H = gca; subplot(2,1,1);
set(H, 'YLimMode', 'manual'); set(H, 'XLimMode', 'manual');
set(H, 'YLim', get(H,'YLim')); set(H, 'XLim', get(H,'XLim'));
xg = []; yg = [];
for i=1:n
    try
        [xi,yi] = ginput(1);
    catch
        return;
    end
    xg = [xg xi]; yg = [yg yi];
    if i == 1
        subplot(2,1,1); hold on;
        plot(H, xg(i),yg(i),'r.'); 
        text(xg(i),yg(i),num2str(i),'FontSize',7);
    end
    if length(xg) > verzögerung
        [signale,avg,dev] = ...
            ThresholdingAlgo(yg,verzögerung,schwelle,einfluss);
        area(xg(verzögerung+1:end),avg(verzögerung+1:end)+schwelle*dev(verzögerung+1:end),...
            'FaceColor',[0.9 0.9 0.9],'EdgeColor','none');
        area(xg(verzögerung+1:end),avg(verzögerung+1:end)-schwelle*dev(verzögerung+1:end),...
            'FaceColor',[1 1 1],'EdgeColor','none');
        plot(xg(verzögerung+1:end),avg(verzögerung+1:end),'LineWidth',1,'Color','cyan');
        plot(xg(verzögerung+1:end),avg(verzögerung+1:end)+schwelle*dev(verzögerung+1:end),...
            'LineWidth',1,'Color','green');
        plot(xg(verzögerung+1:end),avg(verzögerung+1:end)-schwelle*dev(verzögerung+1:end),...
            'LineWidth',1,'Color','green');
        subplot(2,1,2); hold on; title('Signal-Ausgabe');
        stairs(xg(verzögerung+1:end),signale(verzögerung+1:end),'LineWidth',2,'Color','blue');
        ylim([-2 2]); xlim([0 50]); hold off;
    end
    subplot(2,1,1); hold on;
    for j=2:i
        plot(xg([j-1:j]),yg([j-1:j]),'r'); plot(H,xg(j),yg(j),'r.');
        text(xg(j),yg(j),num2str(j),'FontSize',7);
    end
end
end

---

Answer 3

Weitere Algorithmus von Palshikar (2009) gefunden in:

Palshikar, G. (2009). Einfache Algorithmen zur Spitzen­erkennung in Zeitserien. In Proc. 1st Int. Conf. Advanced Data Analysis, Business Analytics and Intelligence (Vol. 122).

Paper kann von hier heruntergeladen werden hier.

Der Algorithmus funktioniert wie folgt:

Algorithmus peak1 // Ein Spitzen­erkennungsalgorithmus, der die Spitzenfunktion S1 verwendet

Eingabe T = x1, x2, …, xN, N // Eingabezeitreihe von N Punkten
Eingabe k // Fenstergröße um den Spitzenwert
Eingabe h // typischerweise 1 <= h <= 3
Ausgabe O // Menge der erkannten Spitzen in T

begin
O = leere Menge // anfangs leer

for (i = 1; i < n; i++) do
// Berechne den Spitzenfunktionswert für jeden der N Punkte in T
a[i] = S1(k,i,xi,T);
end for

Berechne den Mittelwert m' und Standardabweichung s' aller positiven Werte im Array a;

for (i = 1; i < n; i++) do // Entferne lokale Spitzen, die im globalen Kontext "klein" sind
if (a[i] > 0 && (a[i] – m') >( h * s')) then O = O + {xi};
end if
end for

Ordne Spitzen in O in Bezug auf den aufsteigenden Index in T an

// Behalte nur eine Spitze aus jeder Gruppe von Spitzen innerhalb der Entfernung k voneinander bei

für jedes benachbarte Paar von Spitzen xi und xj in O do
if |j – i| <= k then entferne den kleineren Wert von {xi, xj} aus O
end if
end for
end

Vorteile

• Das Papier bietet 5 verschiedene Algorithmen zur Spitzen­erkennung
• Die Algorithmen arbeiten mit den Rohdaten der Zeitreihe (keine Glättung erforderlich)

Nachteile

• Schwierig zu bestimmen k und h im Voraus
• Spitzen können nicht flach sein (wie die dritte Spitze in meinen Testdaten)

Beispiel:

Answer 4

Hier ist eine C-Implementierung von @Jean-Paul's Smoothed Z-score für den Arduino-Mikrocontroller, der verwendet wird, um Beschleunigungsmesswerte zu erfassen und zu entscheiden, ob die Richtung eines Aufpralls von links oder von rechts gekommen ist. Dies funktioniert wirklich gut, da dieses Gerät ein reflektiertes Signal zurückgibt. Hier ist der Eingang zu diesem Peak-Detektionsalgorithmus vom Gerät - der einen Aufprall von rechts und einen Aufprall von links zeigt. Man sieht den initialen Spitzenwert und dann die Schwingung des Sensors.

#include 
#include 
#include 

#define SAMPLE_LENGTH 1000

float stddev(float data[], int len);
float mean(float data[], int len);
void thresholding(float y[], int signals[], int lag, float threshold, float influence);

void thresholding(float y[], int signals[], int lag, float threshold, float influence) {
    memset(signals, 0, sizeof(int) * SAMPLE_LENGTH);
    float filteredY[SAMPLE_LENGTH];
    memcpy(filteredY, y, sizeof(float) * SAMPLE_LENGTH);
    float avgFilter[SAMPLE_LENGTH];
    float stdFilter[SAMPLE_LENGTH];

    avgFilter[lag - 1] = mean(y, lag);
    stdFilter[lag - 1] = stddev(y, lag);

    for (int i = lag; i < SAMPLE_LENGTH; i++) {
        if (fabsf(y[i] - avgFilter[i-1]) > threshold * stdFilter[i-1]) {
            if (y[i] > avgFilter[i-1]) {
                signals[i] = 1;
            } else {
                signals[i] = -1;
            }
            filteredY[i] = influence * y[i] + (1 - influence) * filteredY[i-1];
        } else {
            signals[i] = 0;
        }
        avgFilter[i] = mean(filteredY + i-lag, lag);
        stdFilter[i] = stddev(filteredY + i-lag, lag);
    }
}

float mean(float data[], int len) {
    float sum = 0.0, mean = 0.0;

    int i;
    for(i=0; i

`

Hier ist das Ergebnis mit Einfluss = 0

Nicht großartig, aber hier mit Einfluss = 1

was sehr gut ist.

`

Answer 5

Im Anschluss an @Jean-Pauls vorgeschlagene Lösung habe ich seinen Algorithmus in C# implementiert

public class ZScoreOutput
{
    public List input;
    public List signals;
    public List avgFilter;
    public List filtered_stddev;
}

public static class ZScore
{
    public static ZScoreOutput StartAlgo(List input, int lag, double threshold, double influence)
    {
        // Variablen initialisieren!
        int[] signals = new int[input.Count];
        double[] filteredY = new List(input).ToArray();
        double[] avgFilter = new double[input.Count];
        double[] stdFilter = new double[input.Count];

        var initialWindow = new List(filteredY).Skip(0).Take(lag).ToList();

        avgFilter[lag - 1] = Mean(initialWindow);
        stdFilter[lag - 1] = StdDev(initialWindow);

        for (int i = lag; i < input.Count; i++)
        {
            if (Math.Abs(input[i] - avgFilter[i - 1]) > threshold * stdFilter[i - 1])
            {
                signals[i] = (input[i] > avgFilter[i - 1]) ? 1 : -1;
                filteredY[i] = influence * input[i] + (1 - influence) * filteredY[i - 1];
            }
            else
            {
                signals[i] = 0;
                filteredY[i] = input[i];
            }

            // Rolling-Durchschnitt und Abweichung aktualisieren
            var slidingWindow = new List(filteredY).Skip(i - lag).Take(lag+1).ToList();

            var tmpMean = Mean(slidingWindow);
            var tmpStdDev = StdDev(slidingWindow);

            avgFilter[i] = Mean(slidingWindow);
            stdFilter[i] = StdDev(slidingWindow);
        }

        // In Hilfsklasse kopieren
        var result = new ZScoreOutput();
        result.input = input;
        result.avgFilter       = new List(avgFilter);
        result.signals         = new List(signals);
        result.filtered_stddev = new List(stdFilter);

        return result;
    }

    private static double Mean(List list)
    {
        // Einfache Hilfsfunktion!
        return list.Average();
    }

    private static double StdDev(List values)
    {
        double ret = 0;
        if (values.Count() > 0)
        {
            double avg = values.Average();
            double sum = values.Sum(d => Math.Pow(d - avg, 2));
            ret = Math.Sqrt((sum) / (values.Count() - 1));
        }
        return ret;
    }
}

Beispielverwendung:

var input = new List {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;
double threshold = 5.0;
double influence = 0.0;

var output = ZScore.StartAlgo(input, lag, threshold, influence);