Wie korreliere ich zwei Zeitreihen mit Lücken und unterschiedlichen Zeitbasen?

Ich habe diese Frage bei StackOverflow gestellt und wurde empfohlen, sie hier zu stellen.

Ich habe zwei Zeitreihen von 3D-Beschleunigungsmesserdaten, die unterschiedliche Zeitbasen haben (Uhren wurden zu unterschiedlichen Zeiten gestartet, mit einem sehr geringen Kriechen während der Abtastzeit) sowie viele Lücken unterschiedlicher Größe (aufgrund von Verzögerungen beim Schreiben zum Trennen) Flash-Geräte).

Die Beschleunigungsmesser, die ich benutze, sind die preiswerten GCDC X250-2 . Ich betreibe die Beschleunigungsmesser mit der höchsten Verstärkung, sodass die Daten ein erhebliches Grundrauschen aufweisen.

Die Zeitreihen haben jeweils ungefähr 2 Millionen Datenpunkte (über eine Stunde bei 512 Abtastungen / s) und enthalten ungefähr 500 interessierende Ereignisse, wobei ein typisches Ereignis 100-150 Abtastwerte (jeweils 200-300 ms) umfasst. Viele dieser Ereignisse sind von Datenausfällen während Flash-Schreibvorgängen betroffen.

Die Daten sind also nicht makellos und auch nicht sehr hübsch. Meine Augapfelinspektion zeigt jedoch, dass sie eindeutig die Informationen enthält, an denen ich interessiert bin. (Ich kann bei Bedarf Diagramme veröffentlichen.)

Die Beschleunigungsmesser befinden sich in ähnlichen Umgebungen, sind jedoch nur mäßig gekoppelt, was bedeutet, dass ich anhand des Auges erkennen kann, welche Ereignisse mit den einzelnen Beschleunigungsmessern übereinstimmen. Bisher war dies in der Software jedoch nicht erfolgreich. Aufgrund physikalischer Einschränkungen sind die Geräte auch in unterschiedlichen Ausrichtungen montiert, wobei die Achsen nicht übereinstimmen, aber sie sind so nah an der Orthogonalität, wie ich sie herstellen könnte. So werden beispielsweise für 3-Achsen-Beschleunigungsmesser A und B + Ax auf -By (oben-unten), + Az auf -Bx (links-rechts) und + Ay auf -Bz (vorne-hinten) abgebildet. .

Mein ursprüngliches Ziel ist es, Schockereignisse auf der vertikalen Achse zu korrelieren, obwohl ich schließlich a) die Achsenabbildung automatisch ermitteln, b) die Aktivität auf den abgebildeten Assen korrelieren und c) Verhaltensunterschiede zwischen den beiden Beschleunigungsmessern (z. B. Verdrehen) extrahieren möchte oder biegen).

Die Art der Zeitreihendaten macht Pythons numpy.correlate () unbrauchbar. Ich habe mir auch R's Zoo-Paket angesehen, bin aber damit nicht weitergekommen. Ich habe in verschiedenen Bereichen der Signalanalyse nach Hilfe gesucht, aber keine Fortschritte erzielt.

Hat jemand eine Ahnung, was ich tun kann oder welche Ansätze ich erforschen sollte?

Update 28. Februar 2011: Hier wurden einige Diagramme hinzugefügt , die Beispiele für die Daten zeigen.

time-series correlation unevenly-spaced-time-series BobC
quelle

@ BobC, vielleicht kann einer der Moderatoren Ihren Beitrag auf diese Site migrieren lassen. Das wäre am vernünftigsten. Verwenden Sie für Ihre technischen Fragen zunächst die FFT, um die Korrelation durchzuführen? Dies sollte für 2 Millionen Datenpunkte auf einem halbwegs anständigen Computer möglich sein. Ihr Signal-Rausch-Verhältnis sieht relativ hoch aus, daher sollten Sie im Geschäft sein. Ein schneller und schmutziger Schnitt wäre, die fehlenden Daten entweder mit der letzten verfügbaren Stichprobe oder mit Nullen zu füllen. Das Kriechen aufgrund von Abtastintervallunterschieden kann das schwierigste "Merkmal" Ihrer Daten sein, mit dem Sie sich befassen müssen.

Kardinal

@cardinal: Ich habe in der Tat eine FFT versucht, nur um Müll zu bekommen. Die in den Daten leicht sichtbaren "interessanten" Merkmale sind vom Rauschen in der FFT nicht zu unterscheiden. Ich habe FFTs jedoch nur für den gesamten Datensatz durchgeführt: Möglicherweise würde eine FFT mit beweglichem Fenster bessere Ergebnisse liefern, aber ich konnte noch keinen rechnerisch effizienten Weg finden, um sie zu implementieren. Ich vermute, dass eine Wavelet-Transformation helfen könnte, aber ich bin damit nicht vertraut (lerne aber langsam davon).

BobC

@BobC, was ich damit meinte, haben Sie eine FFT-basierte Implementierung zur Berechnung der Korrelation in Betracht gezogen? Die direkte Faltung ist , aber eine FFT-basierte Implementierung würde dieses reduzieren , was es möglich macht. Wenn Sie sich die FFT selbst ansehen, ist Ihre Frequenzauflösung mit 2 Millionen Datenpunkten sehr hoch. Jedes Abtastkriechen und andere Dinge müssen das Signal pro Frequenz auswaschen. Sie sollten jedoch in der Lage sein, über viele Bins zu aggregieren, um das Signal aus dem Rauschen herauszuholen. So etwas wie ein Welch-Ansatz oder vielleicht eine benutzerdefinierte Fenstertechnik.

O (n^{2})

$O(n^2)$

O (n \log n)

$O(n \log n)$

Kardinal

@BobC, auf den ersten Blick scheint es, als könnte eine Variante eines Algorithmus zum Überlappen und Hinzufügen oder Überlappen und Speichern verwendet werden, um eine Schiebefenster-FFT durchzuführen. Das Verschieben der Samples innerhalb eines Fensters stellt nur eine Phasenverschiebung dar. Sie müssen also nur die Samples kompensieren, die am linken Ende "abfallen" und die am rechten Ende "hereinkommen".

Kardinal

Hallo, ich habe eine ähnliche Frage. Ich habe 2 Zeitreihen, die jeweils durch eine Matrix dargestellt werden, deren erste Spalte den Werten und deren zweite Spalte der Zeitdifferenz (seit dem vorherigen Wert) entspricht. Wie finde ich die Korrelation zwischen diesen beiden Matrizen? Ich habe versucht, xcorr2 () zu machen, aber es scheint nicht richtig zu sein, und xcorr würde wahrscheinlich nur die Korrelation mit den zu berücksichtigenden Werten berechnen, aber ich möchte auch die Zeit berücksichtigen. Ich bin hier wirklich verwirrt. Wird eine FFT helfen? Wie würden Sie vorschlagen, dass ich es mache?

Antworten:

Die Frage betrifft die Berechnung der Korrelation zwischen zwei unregelmäßig abgetasteten Zeitreihen (eindimensionale stochastische Prozesse) und deren Verwendung, um den Zeitversatz zu finden, bei dem sie maximal korreliert sind (ihre "Phasendifferenz").

Dieses Problem wird in der Zeitreihenanalyse normalerweise nicht behandelt, da davon ausgegangen wird, dass Zeitreihendaten systematisch (in regelmäßigen Zeitabständen) erfasst werden. Es ist eher die Provinz der Geostatistik , die sich mit den mehrdimensionalen Verallgemeinerungen von Zeitreihen befasst. Der archetypische geostatistische Datensatz besteht aus Messungen geologischer Proben an unregelmäßig verteilten Orten.

Bei unregelmäßigen Abständen variieren die Abstände zwischen Standortpaaren: Es dürfen keine zwei Abstände gleich sein. Die Geostatistik überwindet dies mit dem empirischen Variogramm . Dies berechnet einen "typischen" (oft den Mittelwert oder Median) Wert von 2/2 - die "Semivarianz" - wobei einen gemessenen Wert am Punkt und der Abstand zwischen und ist darauf beschränkt, innerhalb eines Intervalls zu liegen, das als "Verzögerung" bezeichnet wird. Wenn wir annehmen, dass der Prozess stationär ist und eine Kovarianz aufweist, entspricht die Erwartung der Semivarianz der maximalen Kovarianz (gleich für jedes ) abzüglich der Kovarianz zwischen $(z(p) - z(q))^2 / 2$ $z(p)$ $p$ $p$ $q$ $Z$ $Var(Z(p))$ $p$ $Z(p)$ und . Dieses Binning in Verzögerungen bewältigt das Problem des unregelmäßigen Abstands. $Z(q)$

Wenn an jedem Punkt ein geordnetes Messpaar wird, kann man auf ähnliche Weise das empirische Kreuzvariogramm zwischen den und berechnen und dadurch die Kovarianz bei jeder Verzögerung abschätzen . Sie möchten die eindimensionale Version des Kreuzvariogramms. Die R- Pakete gstat und sgeostat schätzen unter anderem Kreuzvariogramme . Machen Sie sich keine Sorgen, dass Ihre Daten eindimensional sind. Wenn die Software nicht direkt mit ihnen arbeitet, geben Sie einfach eine konstante zweite Koordinate ein, damit sie zweidimensional erscheinen. $(z(p), w(p))$ $z$ $w$

Mit zwei Millionen Punkten sollten Sie kleine Abweichungen von der Stationarität erkennen können. Es ist möglich, dass die Phasendifferenz zwischen den beiden Zeitreihen auch über die Zeit variiert. Bewältigen Sie dies, indem Sie das Kreuzvariogramm für verschiedene Fenster, die über den Zeitraum verteilt sind, separat berechnen.

@cardinal hat die meisten dieser Punkte bereits in Kommentaren angesprochen. Der Hauptbeitrag dieser Antwort besteht darin, auf die Verwendung von Paketen für räumliche Statistiken hinzuweisen, um Ihre Arbeit für Sie zu erledigen, und Techniken der Geostatistik zu verwenden, um diese Daten zu analysieren. Beachten Sie hinsichtlich der Recheneffizienz, dass die vollständige Faltung (Kreuzvariogramm) nicht benötigt wird: Sie benötigen nur ihre Werte in der Nähe der Phasendifferenz. Dies macht den Aufwand , nicht , wobei die Anzahl der zu berechnenden Verzögerungen ist, so dass dies möglicherweise auch mit sofort einsatzbereiter Software möglich ist. Wenn nicht, ist der direkte Faltungsalgorithmus einfach zu implementieren. $O(nk)$ $O(n^2)$ $k$

whuber
quelle

@whuber, gute Kommentare und Vorschläge. Wenn ich die Frage richtig lese, glaube ich, dass ein Hauptanliegen die Unsicherheit des Zeitpunkts der Probenahme ist. Dies kann sich ein wenig von dem typischen geostatistischen Rahmen unterscheiden, bei dem der Abstand meines Erachtens unregelmäßig ist, aber immer noch als bekannt angenommen wird (zumindest mit hoher Präzision). Ich denke, ein grobes Modell ist, wenn der te Punkt der Serie eins zum Zeitpunkt , für festes dann liegt der te Punkt der Serie 2 bei wobei ist wahrscheinlich in der Größenordnung von ein paar Millisekunden und ist wahrscheinlich winzig.

n

$n$

t_{n} = n t

$t_n = n t$

t

$t$

n

$n$

τ_{n} = t_{n} + α + β n

$\tau_n = t_n + \alpha + \beta n$

α

$\alpha$

β

$\beta$

Kardinal

@ Cardinal Ich habe das nicht von der Frage bekommen. Ich kann mir keine Methode zur Schätzung von , die nicht sehr rechenintensiv wäre. Vielleicht durch Aufteilen der Zeitreihen in Gruppen, in denen der Nettoeffekt von vernachlässigbar ist?

β

$\beta$

β

$\beta$

whuber

@whuber, @BobC, ich mache eine halbherzige Vermutung, basierend auf früheren Erfahrungen mit ähnlichen Problemen und Problemen. Die meisten Ansätze, die ich gesehen habe, sind rechenintensiv und beeindrucken nicht. Ein Versuch könnte über etwas wie dynamisches Time Warping oder das, was Ramsay und Silverman Kurvenregistrierung nennen , erfolgen . Ob beides mit diesem Größen-Datensatz machbar wäre, ist mir unklar.

Kardinal

Ich brauche ein bisschen, um mein Gehirn darum zu wickeln. Ich beginne mit den Beispielen in den von Ihnen erwähnten R-Paketen.

BobC

@ BobC, liegt das grobe Modell, das ich für die Timing-Asynchronität angegeben habe, nahe an dem, was Sie haben? Ich denke, es ist "zufälliger Anfangsversatz" + "linearer Fehler", wobei letzterer auf einen kleinen konstanten Unterschied im Abtastintervall zwischen Ihren beiden Geräten zurückzuführen ist. Dann gibt es einen zusätzlichen kleinen zufälligen Fehler, beispielsweise aufgrund der Interrupt-Verarbeitung von zwei verschiedenen uCs.

Kardinal