Znajdowanie punktu zmiany w danych z częściowej funkcji liniowej

Pozdrowienia,

Przeprowadzam badania, które pomogą określić rozmiar obserwowanej przestrzeni i czas, jaki upłynął od Wielkiego Wybuchu. Mam nadzieję, że możesz pomóc!

Mam dane zgodne z częściową funkcją liniową, na której chcę wykonać dwie regresje liniowe. Jest punkt, w którym nachylenie i punkt przecięcia zmieniają się i muszę (napisać program) znaleźć ten punkt.

Myśli?

mcpPakiet może to zrobić. Powiedz, że masz dane

Najpierw symulujmy niektóre dane:

df = data.frame(x = 1:100,
                y = c(rnorm(40, 10 + (1:40)*0.5),
                      rnorm(60, 10 + 40*0.5 -8 + (1:60)*0.2)))

Zobaczmy teraz, czy możemy odzyskać punkt zmiany przy 40 (i wartościach parametrów) za pomocą mcp:

model = list(
  y ~ 1 + x,  # linear segment
  ~ 1 + x  # another linear segment
)
library(mcp)
fit = mcp(model, df)

Działka Szare linie są losowymi losowaniami z dopasowania, co pokazuje, że odzwierciedla trend. Niebieska krzywa jest szacunkową lokalizacją punktu zmiany:

Zobaczmy szacunki poszczególnych parametrów. int_są punktami przecięcia, x_są nachyleniami na x i cp_są punktami zmiany:

summary(fit)

Population-level parameters:
    name  mean lower upper Rhat n.eff
    cp_1 40.48 40.02 41.00    1  2888
   int_1 11.12  9.11 13.17    1   778
   int_2 21.72 20.09 23.49    1   717
 sigma_1  3.23  2.76  3.69    1  5343
     x_1  0.46  0.36  0.54    1   724
     x_2  0.21  0.16  0.26    1   754

Oświadczenie: Jestem deweloperem mcp.

Strucchange pakietu R może ci pomóc. Spójrz na winietę, ma ładny przegląd, jak rozwiązać podobne problemy.

Jeśli liczba punktów nie jest zbyt duża, możesz wypróbować wszystkie możliwości. Załóżmy, że punkty są , gdzie . Następnie możesz zapętlić za pomocą od do i dopasować dwie linie do obu i . Na koniec wybierasz dla którego suma kwadratów reszt dla obu linii jest minimalna.$X_i=(x_i,y_i)$$i=1,..,N$$j$$2$$N-2$$\{X_1,...,X_j\}$$\{X_{(j+1)},...,X_N\}$$j$

Jest to problem z wykrywaniem punktu zmiany (offline). Nasza poprzednia dyskusja zawiera odniesienia do artykułów w czasopismach i kodu R. Najpierw spójrz na „model partycji produktu” Barry'ego i Hartigana , ponieważ obsługuje on zmiany nachylenia i ma wydajne implementacje.

Również pakiet podzielony na segmenty pomógł mi w przeszłości z podobnymi problemami.

Zbudowałem na podstawie odpowiedzi mbq, że szukając wszystkich możliwości. Ponadto robię to:

• Sprawdź znaczenie dwóch modeli częściowych, aby upewnić się, że współczynniki są znaczące
• Sprawdź różnicę do sumy kwadratów reszt dla pełnego modelu
• Potwierdź wizualnie mój model (upewnij się, że to nie jest nonsens)

Po co sprawdzać znaczenie? Wynika to z faktu, że punkt z minimalnym SSE nie ma znaczenia, jeśli któryś z modeli cząstkowych bardzo źle pasuje do danych. Może się to zdarzyć w przypadku dwóch wysoce skorelowanych zmiennych bez wyraźnego punktu przerwania, w którym zmieniają się nachylenia.

Sprawdźmy to proste podejście w prostym przypadku testowym:

x <- c(-50:50)
y <- abs(x)
plot(x,y,pch=19)

Punkt przerwania jest oczywiście zerowy. Użyj następującego skryptu R:

f <- function(x, y)
{
    d <- data.frame(x=x, y=y)
    d <- d[order(x),]
    r <- data.frame(k=rep(0,length(x)-4), sums=rep(0,length(x)-4))

    plm <- function(i)
    {
        d1 <- head(d,i)
        d2 <- tail(d,-i)

        # Make sure we've divided the region perfectly        
        stopifnot(nrow(d1)+nrow(d2) == nrow(d))

        m1 <- lm(y~x, data=d1)
        m2 <- lm(y~x, data=d2)

        r <- list(m1, m2)
        r
    }

    lapply(2:(nrow(d)-3), function(i)
    {
        r$k[i-2] <<- d[i,]$x

        # Fit two piecewise linear models
        m <- plm(i)

        # Add up the sum of squares for residuals
        r$sums[i-2] <<- sum((m[[1]]$residuals)^2) + sum((m[[2]]$residuals)^2)
    })

    b <- r[which.min(r$sums),]    
    b
}

Dopasuj częściowe modele liniowe do wszystkich możliwych kombinacji:

f(x,y)
   k sums
   0    0

Jeśli sprawdzimy współczynniki dla dwóch optymalnych modeli, będą one bardzo znaczące. Ich R2 również będzie bardzo wysoki.