Dane mają dwa trendy; jak wydobywać niezależne linie trendów?

Mam zestaw danych, które nie są uporządkowane w żaden szczególny sposób, ale kiedy są wyraźnie przedstawione, mają dwa wyraźne trendy. Prosta regresja liniowa nie byłaby w tym przypadku wystarczająca ze względu na wyraźne rozróżnienie między dwiema seriami. Czy istnieje prosty sposób na uzyskanie dwóch niezależnych liniowych linii trendu?

Dla przypomnienia korzystam z Pythona i dość dobrze czuję się w programowaniu i analizie danych, w tym w uczeniu maszynowym, ale jestem skłonny przejść do R, jeśli jest to absolutnie konieczne.

Aby rozwiązać problem, dobrym rozwiązaniem jest zdefiniowanie modelu probabilistycznego, który pasuje do założeń dotyczących zestawu danych. W twoim przypadku prawdopodobnie potrzebujesz kombinacji modeli regresji liniowej. Można utworzyć model „mieszaniny regresorów” podobny do modelu mieszanki gaussowskiej, łącząc różne punkty danych z różnymi składnikami mieszanki.

Włączyłem trochę kodu, aby zacząć. Kod implementuje algorytm EM dla mieszanki dwóch regresorów (powinno być względnie łatwe do rozszerzenia na większe mieszaniny). Kod wydaje się być dość niezawodny dla losowych zestawów danych. Jednak w przeciwieństwie do regresji liniowej modele mieszane mają cele niewypukłe, dlatego w przypadku prawdziwego zestawu danych może być konieczne przeprowadzenie kilku prób z różnymi losowymi punktami początkowymi.

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

W innym miejscu tego wątku użytkownik1149913 zapewnia świetne porady (definiuje model probabilistyczny) i koduje potężne podejście (szacowanie EM). Do rozwiązania pozostają dwie kwestie:

1. Jak radzić sobie z odstępstwami od modelu probabilistycznego (które są bardzo widoczne w danych z lat 2011–2012 i nieco widoczne w falowaniu punktów o mniejszym nachyleniu).

2. Jak zidentyfikować dobre wartości początkowe dla algorytmu EM (lub dowolnego innego algorytmu).

Aby rozwiązać problem nr 2, rozważ użycie transformacji Hougha . Jest to algorytm wykrywania cech, który w celu znalezienia liniowych ciągów cech może być skutecznie obliczony jako transformacja Radona .

$x$$y$$x,y$w transformacji Hougha. Kiedy elementy na oryginalnym wykresie spadają wzdłuż wspólnej linii lub blisko niej, wówczas zbiory krzywych, które tworzą w transformacie Hougha, zwykle mają wspólne przecięcie odpowiadające tej wspólnej linii. Znajdując te punkty o największej intensywności w transformacji Hougha, możemy odczytać dobre rozwiązania pierwotnego problemu.

Aby zacząć od tych danych, najpierw wyciąłem elementy pomocnicze (osie, znaczniki i etykiety) i dla pewności wyciąłem oczywiście odległe punkty w prawym dolnym rogu i posypałem wzdłuż dolnej osi. (Gdy te elementy nie są wykadrowane, procedura nadal działa dobrze, ale wykrywa również osie, ramki, liniowe sekwencje kleszczy, liniowe sekwencje etykiet, a nawet punkty leżące sporadycznie na dolnej osi!)

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}]

(To i reszta kodu są w Mathematica .)

Każdej kropce na tym zdjęciu odpowiada wąski zakres krzywych w transformacji Hougha, widoczny tutaj. Są to fale sinusoidalne:

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

To sprawia, że ​​wizualnie uwidacznia się sens, w którym pytanie jest problemem klastrowania linii : transformacja Hougha redukuje go do problemu klastrowania punktowego , do którego możemy zastosować dowolną metodę grupowania.

W tym przypadku grupowanie jest tak jasne, że wystarczy proste przetwarzanie końcowe transformacji Hougha. Aby zidentyfikować lokalizacje o największej intensywności w transformacji, zwiększyłem kontrast i rozmazałem transformację w promieniu około 1%: jest to porównywalne ze średnicą punktów wykresu na oryginalnym obrazie.

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

Próg wyniku zawęził go do dwóch drobnych plamek, których centroidy racjonalnie identyfikują punkty o największej intensywności: szacują dopasowane linie.

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

$0.777$

Lewa strona obrazu odpowiada kierunkowi 0 stopni (poziomo) i, gdy patrzymy od lewej do prawej, kąt ten wzrasta liniowo do 180 stopni. Interpolując, obliczam, że dwie plamy są wyśrodkowane odpowiednio w 19 i 57,1 stopniach. Możemy również odczytać przecięcia z pozycji pionowych obiektów blob. Ta informacja daje początkowe pasowania:

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}

W podobny sposób można obliczyć przecięcia odpowiadające tym zboczom, dając następujące pasowania:

(Czerwona linia odpowiada małej różowej kropce na poprzednim zdjęciu, a niebieska linia odpowiada większej kropli wody).

W dużym stopniu podejście to automatycznie rozwiązało pierwszy problem: odchylenia od liniowości rozmazują punkty o największej intensywności, ale zazwyczaj nie zmieniają ich znacznie. Szczerze mówiąc, odległe punkty przyczynią się do niskiego poziomu hałasu podczas transformacji Hougha, który zniknie podczas procedur przetwarzania końcowego.

W tym momencie można podać te szacunki jako wartości początkowe dla algorytmu EM lub minimalizatora prawdopodobieństwa (który przy dobrych oszacowaniach szybko się zbiegnie). Lepiej byłoby jednak użyć solidnego estymatora regresji, takiego jak iteracyjnie przeważone najmniejsze kwadraty . Jest w stanie zapewnić wagę regresji do każdego punktu. Niskie ciężary wskazują, że punkt nie „należy” do linii. W razie potrzeby wykorzystaj te ciężary, aby przypisać każdy punkt do właściwej linii. Następnie, po sklasyfikowaniu punktów, możesz użyć zwykłych najmniejszych kwadratów (lub dowolnej innej procedury regresji) osobno na dwóch grupach punktów.

Znalazłem to pytanie powiązane z innym pytaniem . Przeprowadziłem akademickie badania tego rodzaju problemu. Proszę sprawdzić moją odpowiedź „Najmniejszy pierwiastek kwadratowy”? Metoda dopasowania z wieloma minimami, aby uzyskać więcej informacji.

Podejście Whubera oparte na transformacji Hougha jest bardzo dobrym rozwiązaniem dla prostych scenariuszy, jak ten, który podałeś. Pracowałem nad scenariuszami z bardziej złożonymi danymi, takimi jak to:

Wspólnie z moimi współautorami określiłem ten problem jako „powiązanie danych”. Kiedy próbujesz go rozwiązać, główny problem jest zazwyczaj kombinatoryczny z powodu wykładniczej ilości możliwych kombinacji danych.

Mamy publikację „ Nakładające się mieszaniny procesów gaussowskich dla problemu powiązania danych ”, w której podeszliśmy do ogólnego problemu krzywych N za pomocą iteracyjnej techniki, dającej bardzo dobre wyniki. Kod Matlaba można znaleźć w artykule.

[Aktualizacja] Implementację techniki OMGP w języku Python można znaleźć w bibliotece GPClust .

Mam inny artykuł, w którym złagodziliśmy problem, aby uzyskać wypukły problem optymalizacji, ale nie został jeszcze zaakceptowany do publikacji. Jest specyficzny dla 2 krzywych, więc idealnie działałby na twoich danych. Daj mi znać jeśli jesteś zainteresowany.

user1149913 ma doskonałą odpowiedź (+1), ale wydaje mi się, że twoje zbieranie danych rozpadło się pod koniec 2011 roku, więc będziesz musiał odciąć tę część swoich danych, a następnie kilka razy uruchomić różne rzeczy z innym losowym współczynniki początkowe, aby zobaczyć, co otrzymujesz.

Jednym prostym sposobem na zrobienie tego jest rozdzielenie danych na dwa zestawy za pomocą oka, a następnie użycie dowolnej techniki modelu liniowego, do której przywykłeś. W R byłaby to lmfunkcja.

Lub dopasuj dwie linie do oka. W R byś to ablinezrobił.

Dane są pomieszane, mają wartości odstające i rozpadają się na końcu, ale naocznie ma dwie dość oczywiste linie, więc nie jestem pewien, czy warto stosować wymyślną metodę.