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

34

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.

wprowadź opis zdjęcia tutaj

jbbiomed
źródło
6
Najlepsza odpowiedź, jaką do tej pory mam, to wydrukować to na papierze
milimetrowym
Może uda ci się obliczyć zbocza parami i zgrupować je w dwie „grupy zboczy”. Jednak to się nie powiedzie, jeśli masz dwa równoległe trendy.
Thomas Jungblut,
1
Nie mam z tym osobistego doświadczenia, ale myślę, że warto sprawdzić statsmodels . Statystycznie regresja liniowa z interakcją dla grupy byłaby wystarczająca (chyba że mówisz, że masz niezgrupowane dane, w którym to przypadku jest to trochę bardziej włochate ...)
Matt Parker
1
Niestety nie są to dane dotyczące efektów, ale dane dotyczące użytkowania, a także wyraźnie wykorzystanie z dwóch oddzielnych systemów zmieszanych w tym samym zestawie danych. Chcę być w stanie opisać dwa wzorce użytkowania, ale nie mogę wrócić i zebrać danych, ponieważ stanowią one informacje zebrane przez klienta o wartości około 6 lat.
jbbiomed
2
Dla pewności: twój klient nie ma żadnych dodatkowych danych, które wskazywałyby, które pomiary pochodzą z jakiej populacji? Jest to 100% danych, które Ty lub Twój klient macie lub możecie znaleźć. Również w 2012 r. Wygląda na to, że Twój zbiór danych się rozpadł lub jeden lub oba systemy upadły na podłogę. Zastanawiam się, czy trendy do tego momentu mają duże znaczenie.
Wayne

Odpowiedzi:

30

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()
użytkownik1149913
źródło
25

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 .

xyx,yw 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 .)

Przycięty obraz

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

Przekształcenie Hougha

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

Niewyraźna transformacja

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

Progowa transformacja binarna

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:

Dopasowane linie

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

Whuber
źródło
1
Zdjęcia mówią tysiąc słów, a ty masz 5. To niesamowita praca z szybkiego wykresu, który stworzyłem tylko na potrzeby tego pytania! Sława!
jbbiomed
2
Przekształcenie Hougha jest szeroko stosowane w polu widzenia komputerowego do identyfikacji linii prostych na obrazie. Dlaczego nie powinien być stosowany również w statystykach? ;)
Lucas Reis
xy
Tak. Wyobraź sobie na przykład liczbę wartości odstających zaangażowanych w porównywanie dwóch obrazów w celu wykrycia, czy pochodzą one od tego samego obiektu. A przede wszystkim wyobraź sobie, że musisz to robić w czasie rzeczywistym. „Szybkość” jest bardzo ważnym czynnikiem w wizji komputerowej, a nie tak ważnym w statystyce.
Lucas Reis,
@RoyalTS Dziękujemy za zwrócenie uwagi na potrzebę naprawy jednego z fragmentów kodu. Zanim znalazłem twoją sugerowaną zmianę, została ona odrzucona (poprawnie, ponieważ nie była w porządku, ale nieważne, że: jestem wdzięczny, że zauważyłeś błąd). Naprawiłem to, usuwając odniesienie do rotation, które pierwotnie było zerowane i dlatego nie robiło różnicy.
whuber
15

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:

problem powiązania danych - zestaw danych cukierków

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.

Steven
źródło
1
Z przykrością widzę, że przez ponad dwa lata nikt inny nie ocenił tej oryginalnej i wartościowej odpowiedzi. Czy w międzyczasie zaakceptowano ostatni wspomniany artykuł?
whuber
1
Artykuł rzeczywiście został przyjęty zaledwie kilka miesięcy temu. Możesz go pobrać tutaj gtas.unican.es/pub/378 . To właściwie dość rzadki problem (co może tłumaczyć jego brak popularności), ale wciąż udało nam się znaleźć kilka interesujących aplikacji. Jeśli chcesz, spójrz na eksperymenty na końcu artykułu.
Steven
2

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ę.

Wayne
źródło