Daten weisen zwei Trends auf; Wie extrahiere ich unabhängige Trendlinien?

34

Ich habe eine Reihe von Daten, die nicht in einer bestimmten Reihenfolge angeordnet sind, aber bei einer klaren Darstellung zwei unterschiedliche Trends aufweisen. Eine einfache lineare Regression wäre hier aufgrund der eindeutigen Unterscheidung der beiden Reihen nicht ausreichend. Gibt es eine einfache Möglichkeit, die beiden unabhängigen linearen Trendlinien zu ermitteln?

Ich benutze Python und bin einigermaßen vertraut mit Programmierung und Datenanalyse, einschließlich maschinellem Lernen, bin aber bereit, bei Bedarf auf R umzusteigen.

Bildbeschreibung hier eingeben

jbbiomed
quelle
6
Die beste Antwort, die ich bisher habe, ist, sie auf
Millimeterpapier
Vielleicht können Sie paarweise Steigungen berechnen und zu zwei "Steigungsclustern" zusammenfassen. Dies schlägt jedoch fehl, wenn Sie zwei parallele Trends haben.
Thomas Jungblut
1
Ich habe keine persönlichen Erfahrungen damit, aber ich denke, dass es sich lohnen würde, die Statistikmodelle zu testen . Statistisch gesehen wäre eine lineare Regression mit einer Interaktion für die Gruppe angemessen (es sei denn, Sie haben nicht gruppierte Daten, in diesem Fall ist das etwas haariger ...)
Matt Parker
1
Leider handelt es sich hierbei nicht um Effektdaten, sondern um Nutzungsdaten, und die Nutzung von zwei separaten Systemen wird eindeutig in denselben Datensatz gemischt. Ich möchte in der Lage sein, die beiden Verwendungsmuster zu beschreiben, kann mich aber nicht an Daten erinnern, da dies einen Wert von 6 Jahren darstellt, der von einem Kunden gesammelt wurde.
jbbiomed
2
Nur um sicherzugehen: Ihr Kunde hat keine zusätzlichen Daten, aus denen hervorgeht, welche Messungen aus welcher Grundgesamtheit stammen. Dies sind 100% der Daten, die Sie oder Ihr Kunde haben oder finden können. 2012 scheint entweder Ihre Datenerfassung auseinandergefallen zu sein oder eines oder beide Ihrer Systeme sind durch den Boden gefallen. Ich frage mich, ob Trendlinien bis zu diesem Punkt eine große Rolle spielen.
Wayne

Antworten:

30

Ein guter Ansatz zur Lösung Ihres Problems besteht darin, ein Wahrscheinlichkeitsmodell zu definieren, das mit den Annahmen zu Ihrem Dataset übereinstimmt. In Ihrem Fall möchten Sie wahrscheinlich eine Mischung aus linearen Regressionsmodellen. Sie können ein "Gemisch von Regressoren" -Modell erstellen, das einem Gaußschen Gemischmodell ähnelt, indem Sie verschiedene Datenpunkte mit verschiedenen Gemischkomponenten verknüpfen.

Ich habe Code eingefügt, um Ihnen den Einstieg zu erleichtern. Der Code implementiert einen EM-Algorithmus für eine Mischung aus zwei Regressoren (eine Erweiterung auf größere Mischungen sollte relativ einfach sein). Der Code scheint für zufällige Datensätze ziemlich robust zu sein. Im Gegensatz zur linearen Regression haben Mischungsmodelle jedoch nicht konvexe Ziele. Daher müssen Sie für einen realen Datensatz möglicherweise einige Versuche mit verschiedenen zufälligen Startpunkten durchführen.

import numpy as np
import matplotlib.pyplot as plt 
import scipy.linalg as lin

#generate some random data
N=100
x=np.random.rand(N,2)
x[:,1]=1

w=np.random.rand(2,2)
y=np.zeros(N)

n=int(np.random.rand()*N)
y[:n]=np.dot(x[:n,:],w[0,:])+np.random.normal(size=n)*.01
y[n:]=np.dot(x[n:,:],w[1,:])+np.random.normal(size=N-n)*.01


rx=np.ones( (100,2) )
r=np.arange(0,1,.01)
rx[:,0]=r

#plot the random dataset
plt.plot(x[:,0],y,'.b')
plt.plot(r,np.dot(rx,w[0,:]),':k',linewidth=2)
plt.plot(r,np.dot(rx,w[1,:]),':k',linewidth=2)

# regularization parameter for the regression weights
lam=.01

def em():
    # mixture weights
    rpi=np.zeros( (2) )+.5

    # expected mixture weights for each data point
    pi=np.zeros( (len(x),2) )+.5

    #the regression weights
    w1=np.random.rand(2)
    w2=np.random.rand(2)

    #precision term for the probability of the data under the regression function 
    eta=100

    for _ in xrange(100):
        if 0:
            plt.plot(r,np.dot(rx,w1),'-r',alpha=.5)
            plt.plot(r,np.dot(rx,w2),'-g',alpha=.5)

        #compute lhood for each data point
        err1=y-np.dot(x,w1)
        err2=y-np.dot(x,w2)
        prbs=np.zeros( (len(y),2) )
        prbs[:,0]=-.5*eta*err1**2
        prbs[:,1]=-.5*eta*err2**2

        #compute expected mixture weights
        pi=np.tile(rpi,(len(x),1))*np.exp(prbs)
        pi/=np.tile(np.sum(pi,1),(2,1)).T

        #max with respect to the mixture probabilities
        rpi=np.sum(pi,0)
        rpi/=np.sum(rpi)

        #max with respect to the regression weights
        pi1x=np.tile(pi[:,0],(2,1)).T*x
        xp1=np.dot(pi1x.T,x)+np.eye(2)*lam/eta
        yp1=np.dot(pi1x.T,y)
        w1=lin.solve(xp1,yp1)

        pi2x=np.tile(pi[:,1],(2,1)).T*x
        xp2=np.dot(pi2x.T,x)+np.eye(2)*lam/eta
        yp2=np.dot(pi[:,1]*y,x)
        w2=lin.solve(xp2,yp2)

        #max wrt the precision term
        eta=np.sum(pi)/np.sum(-prbs/eta*pi)

        #objective function - unstable as the pi's become concentrated on a single component
        obj=np.sum(prbs*pi)-np.sum(pi[pi>1e-50]*np.log(pi[pi>1e-50]))+np.sum(pi*np.log(np.tile(rpi,(len(x),1))))+np.log(eta)*np.sum(pi)
        print obj,eta,rpi,w1,w2

        try:
            if np.isnan(obj): break
            if np.abs(obj-oldobj)<1e-2: break
        except:
            pass

        oldobj=obj

    return w1,w2


#run the em algorithm and plot the solution
rw1,rw2=em()
plt.plot(r,np.dot(rx,rw1),'-r')
plt.plot(r,np.dot(rx,rw2),'-g')

plt.show()
user1149913
quelle
25

An anderer Stelle in diesem Thread bietet user1149913 hervorragende Ratschläge (Definieren eines probabilistischen Modells) und Code für einen leistungsfähigen Ansatz (EM-Schätzung). Zwei Fragen müssen noch geklärt werden:

  1. Wie gehe ich mit Abweichungen vom Wahrscheinlichkeitsmodell um (die in den Daten für 2011-2012 sehr deutlich und in den Unebenheiten der weniger geneigten Punkte etwas deutlich erkennbar sind)?

  2. Identifizieren guter Startwerte für den EM-Algorithmus (oder einen anderen Algorithmus).

Verwenden Sie eine Hough-Transformation, um Nummer 2 zu adressieren . Dies ist ein Merkmalserkennungsalgorithmus, der zum Auffinden linearer Merkmalsabschnitte effizient als Radon-Transformation berechnet werden kann .

xyx,yin der Hough-Transformation. Wenn Merkmale in der ursprünglichen Zeichnung entlang einer gemeinsamen Linie oder nahe genug bei eins liegen, haben die Sammlungen von Kurven, die sie in der Hough-Transformation erzeugen, tendenziell einen gemeinsamen Schnittpunkt, der dieser gemeinsamen Linie entspricht. Indem wir diese Punkte mit der größten Intensität in der Hough-Transformation finden, können wir gute Lösungen für das ursprüngliche Problem ablesen.

Um mit diesen Daten zu beginnen, habe ich zuerst die Hilfsmittel (Achsen, Häkchen und Beschriftungen) und zum guten Teil die offensichtlich äußeren Punkte unten rechts ausgeschnitten und entlang der unteren Achse gestreut. (Wenn das Zeug nicht ausgeschnitten ist, funktioniert die Prozedur immer noch gut, aber es erkennt auch die Achsen, die Frames, die linearen Sequenzen von Ticks, die linearen Sequenzen von Labels und sogar die Punkte, die sporadisch auf der unteren Achse liegen!)

img = Import["http://i.stack.imgur.com/SkEm3.png"]
i = ColorNegate[Binarize[img]]
crop2 = ImageCrop[ImageCrop[i, {694, 531}, {Left, Bottom}], {565, 467}, {Right, Top}]

(Dies und der Rest des Codes befinden sich in Mathematica .)

Zugeschnittenes Bild

Jedem Punkt in diesem Bild entspricht ein enger Kurvenbereich in der Hough-Transformation, der hier sichtbar ist. Sie sind Sinuswellen:

hough2 = Radon[crop2, Method -> "Hough"]  // ImageAdjust

Hough verwandeln

Dies verdeutlicht visuell den Sinn, in dem die Frage ein Linienclustering- Problem ist: Die Hough-Transformation reduziert sie auf ein Punktclustering- Problem, auf das wir jede beliebige Clustering-Methode anwenden können.

In diesem Fall ist die Clusterbildung so eindeutig, dass eine einfache Nachbearbeitung der Hough-Transformation ausreicht. Um Stellen mit der größten Intensität in der Transformation zu identifizieren, habe ich den Kontrast erhöht und die Transformation über einen Radius von etwa 1% verwischt: Das ist vergleichbar mit den Durchmessern der Diagrammpunkte im Originalbild.

blur = ImageAdjust[Blur[ImageAdjust[hough2, {1, 0}], 8]]

Verschwommene Transformation

Durch Schwellenwertbildung wurde das Ergebnis auf zwei winzige Flecken eingegrenzt, deren Zentroide die Punkte mit der größten Intensität identifizieren: Diese schätzen die angepassten Linien.

comp = MorphologicalComponents[blur, 0.777]) // Colorize

0,777

Binäre Transformation mit Schwellenwert

Die linke Seite des Bildes entspricht einer Richtung von 0 Grad (horizontal), und wenn wir von links nach rechts schauen, vergrößert sich dieser Winkel linear auf 180 Grad. Wenn ich interpoliere, berechne ich, dass die beiden Blobs bei 19 bzw. 57,1 Grad zentriert sind. Wir können die Abschnitte auch an den vertikalen Positionen der Blobs ablesen. Diese Information liefert die Anfangsanpassungen:

width = ImageDimensions[blur][[1]];
slopes =  Module[{x, y, z}, ComponentMeasurements[comp, "Centroid"] /. 
          Rule[x_, {y_, z_}] :>  Round[((y - 1/2)/(width - 1))  180., 0.1]
  ]

{19., 57.1}

Auf ähnliche Weise kann man die Abschnitte berechnen, die diesen Steigungen entsprechen, und diese Anpassungen ergeben:

Eingepasste Linien

(Die rote Linie entspricht dem winzigen rosa Punkt im vorherigen Bild und die blaue Linie entspricht dem größeren Aqua-Blob.)

Dieser Ansatz hat sich weitgehend automatisch mit dem ersten Problem befasst: Abweichungen von der Linearität verwischen die Punkte mit der größten Intensität, verschieben sie jedoch in der Regel nicht stark. Ehrlich gesagt tragen äußere Punkte während der gesamten Hough-Transformation zu einem niedrigen Rauschpegel bei, der während der Nachbearbeitungsvorgänge verschwindet.

An diesem Punkt kann man diese Schätzungen als Startwerte für den EM-Algorithmus oder für einen Likelihood-Minimierer (der bei guten Schätzungen schnell konvergiert) bereitstellen. Besser wäre es jedoch, einen robusten Regressionsschätzer wie iterativ gewichtete kleinste Quadrate zu verwenden . Es ist in der Lage , jedem Punkt ein Regressionsgewicht zuzuweisen . Niedrige Gewichte bedeuten, dass ein Punkt nicht zu einer Linie gehört. Nutzen Sie diese Gewichte, wenn Sie möchten, um jedem Punkt die richtige Linie zuzuweisen. Nachdem Sie die Punkte klassifiziert haben, können Sie gewöhnliche kleinste Quadrate (oder ein anderes Regressionsverfahren) für die beiden Punktgruppen separat verwenden.

whuber
quelle
1
Bilder sagen mehr als tausend Worte und Sie haben 5. Dies ist eine unglaubliche Arbeit von einem schnellen Diagramm, das ich nur für den Zweck dieser Frage gemacht habe! Ein dickes Lob!
jbbiomed
2
Hough-Transformation wird im Bereich Computer Vision häufig verwendet, um gerade Linien in einem Bild zu identifizieren. Warum sollte es nicht auch in der Statistik verwendet werden? ;)
Lucas Reis
xy
Ja. Stellen Sie sich zum Beispiel vor, wie viele Ausreißer beim Vergleich zweier Bilder beteiligt sind, um festzustellen, ob sie zum selben Thema gehören. Und vor allem stellen Sie sich vor, Sie müssten es in Echtzeit tun. "Geschwindigkeit" ist ein sehr wichtiger Faktor in der Bildverarbeitung und nicht so wichtig in der Statistik.
Lucas Reis
@ RoyalTS Vielen Dank, dass Sie auf die Notwendigkeit einer Korrektur für einen der Codeausschnitte hingewiesen haben. Als ich Ihre vorgeschlagene Änderung fand, wurde sie abgelehnt (richtig, weil sie nicht ganz richtig war, aber egal: Ich bin dankbar, dass Sie einen Fehler festgestellt haben). Ich habe es durch Entfernen des Verweises auf behoben, der rotationursprünglich auf Null gesetzt war und daher keinen Unterschied machte.
Whuber
15

Ich habe diese Frage mit einer anderen Frage verknüpft . Ich habe tatsächlich akademisch über diese Art von Problem geforscht. Bitte überprüfen Sie meine Antwort "Least Square Root" passend? Eine Anpassungsmethode mit mehreren Minima für weitere Details.

Der auf Hough-Transformation basierende Ansatz von whuber ist eine sehr gute Lösung für einfache Szenarien, wie Sie sie angegeben haben. Ich habe an Szenarien mit komplexeren Daten gearbeitet, wie zum Beispiel:

Datenassoziationsproblem - Süßigkeitendatensatz

Meine Co-Autoren und ich bezeichneten dies als ein "Datenassoziations" -Problem. Wenn Sie versuchen, es zu lösen, ist das Hauptproblem aufgrund der exponentiellen Menge möglicher Datenkombinationen in der Regel kombinatorisch.

Wir haben eine Veröffentlichung " Überlappende Mischungen von Gaußschen Prozessen für das Datenassoziationsproblem " veröffentlicht, in der wir uns dem allgemeinen Problem der N-Kurven mit einer iterativen Technik näherten und sehr gute Ergebnisse erzielten. Sie finden den Matlab-Code im Artikel.

[Update] Eine Python-Implementierung der OMGP-Technik finden Sie in der GPClust-Bibliothek .

Ich habe ein anderes Papier, in dem wir das Problem gelockert haben, um ein konvexes Optimierungsproblem zu erhalten, aber es wurde noch nicht zur Veröffentlichung angenommen. Es ist spezifisch für 2 Kurven, sodass es perfekt für Ihre Daten funktioniert. Lass es mich wissen wenn du interessiert bist.

Steven
quelle
1
Ich bin traurig zu sehen, dass über zwei Jahre niemand diese originelle und wertvolle Antwort gebilligt hat. Wurde in der Zwischenzeit das letzte von Ihnen erwähnte Papier angenommen?
whuber
1
Das Papier wurde tatsächlich erst vor einigen Monaten angenommen. Sie können es hier herunterladen: gtas.unican.es/pub/378 . Dies ist eigentlich ein recht seltenes Problem (was möglicherweise die mangelnde Beliebtheit erklärt), aber wir haben es trotzdem geschafft, einige interessante Anwendungen zu finden. Schauen Sie sich die Experimente am Ende des Artikels an, wenn Sie möchten.
Steven
2

user1149913 hat eine exzellente Antwort (+1), aber es sieht so aus, als ob Ihre Datenerfassung Ende 2011 auseinanderfiel. Sie müssten also diesen Teil Ihrer Daten abschneiden und die Dinge dann immer noch ein paar Mal mit einem anderen Zufallsprinzip ausführen Startkoeffizienten, um zu sehen, was Sie bekommen.

Eine einfache Möglichkeit besteht darin, Ihre Daten per Auge in zwei Gruppen zu unterteilen und dann die gewohnte lineare Modelltechnik zu verwenden. In R wäre es die lmFunktion.

Oder passen Sie zwei Linien mit dem Auge an. In R würden Sie dies verwenden abline.

Die Daten sind durcheinander, haben Ausreißer und fallen am Ende auseinander, aber By-Eye hat zwei ziemlich offensichtliche Linien, daher bin ich mir nicht sicher, ob sich eine ausgefallene Methode lohnt.

Wayne
quelle