Jak wykonać częściową regresję liniową z wieloma nieznanymi węzłami?

Czy są jakieś pakiety do częściowej regresji liniowej, które mogą automatycznie wykrywać wiele węzłów? Dzięki. Kiedy korzystam z pakietu strucchange. Nie mogłem wykryć punktów zmiany. Nie mam pojęcia, jak wykrywa punkty zmiany. Na podstawie wykresów widziałem, że jest kilka punktów, które mogą pomóc mi je wybrać. Czy ktoś mógłby podać tutaj przykład?

Czy MARS miałby zastosowanie? R ma pakiet, earthktóry go implementuje.

Ogólnie rzecz biorąc, trochę dziwne jest dopasowanie czegoś jako liniowego kawałka. Jednak jeśli naprawdę chcesz to zrobić, algorytm MARS jest najbardziej bezpośredni. Będzie budować funkcję jeden węzeł na raz; a następnie zwykle przycina liczbę węzłów, aby zwalczyć nadmiernie ala drzewa decyzyjne. Możesz uzyskać dostęp do algorytmu MARS w R za pośrednictwem earthlub mda. Ogólnie rzecz biorąc, jest on zgodny z GCV, który nie jest tak daleko od pozostałych kryteriów informacyjnych (AIC, BIC itp.)

MARS tak naprawdę nie zapewni ci „optymalnego” dopasowania, ponieważ węzły są hodowane pojedynczo. Naprawdę trudno byłoby dopasować naprawdę „optymalną” liczbę węzłów, ponieważ możliwe kombinacje umieszczania węzłów szybko by wybuchły.

Ogólnie dlatego ludzie zwracają się w kierunku wygładzania splajnów. Większość wygładzających splajnów ma kształt sześcienny, dzięki czemu można oszukać ludzkie oko, by nie zauważył nieciągłości. Byłoby jednak możliwe wykonanie liniowego wygładzania splajnu. Dużą zaletą wygładzania splajnów jest ich pojedynczy parametr do optymalizacji. To pozwala szybko osiągnąć naprawdę „optymalne” rozwiązanie bez konieczności przeszukiwania mnóstwa permutacji. Jeśli jednak naprawdę chcesz szukać punktów przegięcia i masz wystarczająco dużo danych, aby to zrobić, prawdopodobnie najlepszym rozwiązaniem byłby MARS.

Oto przykładowy kod karanych wygładzaniem liniowym splajnów w R:

require(mgcv);data(iris);
gam.test <- gam(Sepal.Length ~ s(Petal.Width,k=6,bs='ps',m=0),data=iris)
summary(gam.test);plot(gam.test);

Rzeczywiste wybrane węzły niekoniecznie korelują z żadnym prawdziwym punktem przegięcia.

Kilka lat temu zaprogramowałem to od zera i mam plik Matlab do wykonywania regresji liniowej na komputerze. Około 1 do 4 punktów przerwania jest obliczeniowo możliwe dla około 20 punktów pomiaru. 5 lub 7 punktów przerwania zaczyna być naprawdę za dużo.

Podejście czysto matematyczne, jak widzę, polega na wypróbowaniu wszystkich możliwych kombinacji, jak sugeruje użytkownik mbq w pytaniu powiązanym z komentarzem pod twoim pytaniem.

Ponieważ wszystkie dopasowane linie są kolejne i przylegają do siebie (bez nakładania się), kombinatoryka podąży za trójkątem Paskal. Gdyby segmenty linii nakładały się między wykorzystanymi punktami danych, uważam, że kombinatoryka podążyłaby za liczbami Stirlinga drugiego rodzaju.

Moim zdaniem najlepszym rozwiązaniem jest wybór kombinacji dopasowanych linii, która ma najniższe odchylenie standardowe wartości korelacji R ^ 2 dopasowanych linii. Spróbuję wyjaśnić na przykładzie. Pamiętaj jednak, że pytanie, ile punktów przerwania należy znaleźć w danych, jest podobne do pytania „Jak długie jest wybrzeże Wielkiej Brytanii?” jak w jednym z artykułów Benoita Mandelbrotsa (matematyka) na temat fraktali. I istnieje kompromis między liczbą punktów przerwania a głębokością regresji.

Teraz do przykładu.

Załóżmy, że mamy idealne dane w funkcji ( i są liczbami całkowitymi):x x $y$$x$$x$$y$

$$\begin{array}{|c|c|c|c|c|c|} \hline &x &y &R^2 line 1 &R^2 line 2 &sum of R^2 values &standard deviation of R^2 \\ \hline &1 &1 &1,000 &0,0400 &1,0400 &0,6788 \\ \hline &2 &2 &1,000 &0,0118 &1,0118 &0,6987 \\ \hline &3 &3 &1,000 &0,0004 &1,0004 &0,7067 \\ \hline &4 &4 &1,000 &0,0031 &1,0031 &0,7048 \\ \hline &5 &5 &1,000 &0,0135 &1,0135 &0,6974 \\ \hline &6 &6 &1,000 &0,0238 &1,0238 &0,6902 \\ \hline &7 &7 &1,000 &0,0277 &1,0277 &0,6874 \\ \hline &8 &8 &1,000 &0,0222 &1,0222 &0,6913 \\ \hline &9 &9 &1,000 &0,0093 &1,0093 &0,7004 \\ \hline &10 &10 &1,000 &-1,978 &1,000 &0,7071 \\ \hline &11 &9 &0,9709 &0,0271 &0,9980 &0,6673 \\ \hline &12 &8 &0,8951 &0,1139 &1,0090 &0,5523 \\ \hline &13 &7 &0,7734 &0,2558 &1,0292 &0,3659 \\ \hline &14 &6 &0,6134 &0,4321 &1,0455 &0,1281 \\ \hline &15 &5 &0,4321 &0,6134 &1,0455 &0,1282 \\ \hline &16 &4 &0,2558 &0,7733 &1,0291 &0,3659 \\ \hline &17 &3 &0,1139 &0,8951 &1,0090 &0,5523 \\ \hline &18 &2 &0,0272 &0,9708 &0,9980 &0,6672 \\ \hline &19 &1 &0 &1,000 &1,000 &0,7071 \\ \hline &20 &2 &0,0094 &1,000 &1,0094 &0,7004 \\ \hline &21 &3 &0,0222 &1,000 &1,0222 &0,6914 \\ \hline &22 &4 &0,0278 &1,000 &1,0278 &0,6874 \\ \hline &23 &5 &0,0239 &1,000 &1,0239 &0,6902 \\ \hline &24 &6 &0,0136 &1,000 &1,0136 &0,6974 \\ \hline &25 &7 &0,0032 &1,000 &1,0032 &0,7048 \\ \hline &26 &8 &0,0004 &1,000 &1,0004 &0,7068 \\ \hline &27 &9 &0,0118 &1,000 &1,0118 &0,6987 \\ \hline &28 &10 &0,04 &1,000 &1,04 &0,6788 \\ \hline \end{array}$$

Te wartości y mają wykres:

Który ma wyraźnie dwa punkty przerwania. Dla celów argumentu obliczymy wartości korelacji R ^ 2 (z formułami komórek Excela (europejski styl przecinek-kropka)):

=INDEX(LINEST(B1:$B$1;A1:$A$1;TRUE;TRUE);3;1)
=INDEX(LINEST(B1:$B$28;A1:$A$28;TRUE;TRUE);3;1)

dla wszystkich możliwych nie nakładających się kombinacji dwóch dopasowanych linii. Wszystkie możliwe pary wartości R ^ 2 mają wykres:

Pytanie brzmi, którą parę wartości R ^ 2 powinniśmy wybrać i jak uogólnić do wielu punktów przerwania, jak podano w tytule? Jednym wyborem jest wybór kombinacji, dla której suma korelacji R-kwadrat jest najwyższa. Kreśląc to, otrzymujemy górną niebieską krzywą poniżej:

Niebieska krzywa, suma wartości R-kwadrat, jest najwyższa pośrodku. Jest to lepiej widoczne z tabeli o wartości jako najwyższej wartości. Jednak moim zdaniem minimum czerwonej krzywej jest dokładniejsze. To znaczy, minimalne odchylenie standardowe wartości R ^ 2 dopasowanych linii regresji powinno być najlepszym wyborem.$1,0455$

Regresja liniowa według kawałków - Matlab - wiele punktów przerwania

Istnieje całkiem niezły algorytm opisany w Tomé i Miranda (1984) .

Proponowana metodologia wykorzystuje podejście najmniejszych kwadratów do obliczenia najlepszego ciągłego zestawu linii prostych, które pasują do danego szeregu czasowego, z zastrzeżeniem szeregu ograniczeń dotyczących minimalnej odległości między punktami przerwania i minimalnej zmiany trendu w każdym punkcie przerwania.

Kod i GUI są dostępne zarówno w Fortran, jak i IDL na ich stronie internetowej: http://www.dfisica.ubi.pt/~artome/linearstep.html

... przede wszystkim musisz to zrobić iteracyjnie i pod pewnymi kryteriami informacyjnymi, takimi jak AIC AICc BIC Cp; ponieważ można uzyskać „idealne” dopasowanie, jeśli liczba węzłów K = liczba punktów danych N, ok. ... najpierw wstaw K = 0; oszacuj L = K + 1 regresji, oblicz na przykład AICc; następnie załóż minimalną liczbę punktów danych w oddzielnym segmencie, powiedzmy L = 3 lub L = 4, ok ... wstaw K = 1; zacznij od L-tych danych jako pierwszego węzła, oblicz SS lub MLE, ... i krok po kroku następny punkt danych jako węzeł, SS lub MLE, aż do ostatniego węzła na danych N-L; wybierz układ z najlepszym dopasowaniem (SS lub MLE) oblicz AICc ... ... wstaw K = 2; ... użyj wszystkich poprzednich regresji (to jest ich SS lub MLE), ale krok po kroku podziel jeden segment na wszystkie możliwe części ... wybierz układ z najlepszym dopasowaniem (SS lub MLE) oblicz AICc ... jeśli ostatni AICc występuje wyżej niż poprzedni: zatrzymaj iteracje! To optymalne rozwiązanie zgodnie z kryterium AICc, ok

Kiedyś natknąłem się na program o nazwie Joinpoint . Na swojej stronie internetowej mówią, że pasuje to do modelu punktów łączenia, w którym „kilka różnych linii jest połączonych razem w„ punktach połączenia ”. I dalej: „Użytkownik podaje minimalną i maksymalną liczbę punktów łączenia. Program rozpoczyna się od minimalnej liczby punktów łączenia (np. 0 punktów łączenia, co jest linią prostą) i sprawdza, czy więcej punktów łączenia jest statystycznie znaczących i należy je dodać do modelu (do tej maksymalnej liczby). ”

NCI używa go do modelowania trendów wskaźników zachorowań na raka, być może pasuje to również do twoich potrzeb.

Aby dopasować do danych, funkcja kawałkowa:

$a_1 , a_2 , p_1 , q_1, p_2 , q_2 , p_3 , q_3$

Na przykład, przy dokładnych danych dostarczonych przez Mats Granvik, wynik jest następujący:

Bez rozproszonych danych ten przykład nie jest bardzo znaczący. Inne przykłady z rozproszonymi danymi pokazano w odnośniku.

Możesz skorzystać z mcppakietu, jeśli znasz liczbę punktów zmiany do wnioskowania. Daje to dużą elastyczność modelowania i wiele informacji o punktach zmiany i parametrach regresji, ale kosztem szybkości.

Witryna mcp zawiera wiele zastosowanych przykładów, np.

library(mcp)

# Define the model
model = list(
  response ~ 1,  # plateau (int_1)
  ~ 0 + time,    # joined slope (time_2) at cp_1
  ~ 1 + time     # disjoined slope (int_3, time_3) at cp_2
)

# Fit it. The `ex_demo` dataset is included in mcp
fit = mcp(model, data = ex_demo)

Następnie możesz wizualizować:

plot(fit)

Lub podsumuj:

summary(fit)

Family: gaussian(link = 'identity')
Iterations: 9000 from 3 chains.
Segments:
  1: response ~ 1
  2: response ~ 1 ~ 0 + time
  3: response ~ 1 ~ 1 + time

Population-level parameters:
    name match  sim  mean lower  upper Rhat n.eff
    cp_1    OK 30.0 30.27 23.19 38.760    1   384
    cp_2    OK 70.0 69.78 69.27 70.238    1  5792
   int_1    OK 10.0 10.26  8.82 11.768    1  1480
   int_3    OK  0.0  0.44 -2.49  3.428    1   810
 sigma_1    OK  4.0  4.01  3.43  4.591    1  3852
  time_2    OK  0.5  0.53  0.40  0.662    1   437
  time_3    OK -0.2 -0.22 -0.38 -0.035    1   834

Oświadczenie: Jestem programistą mcp.