Ermitteln des Änderungspunkts in Daten aus einer stückweise linearen Funktion

Schöne Grüße,

Ich forsche, um die Größe des beobachteten Raums und die seit dem Urknall verstrichene Zeit zu bestimmen. Hoffentlich kannst du helfen!

Ich habe Daten, die einer stückweise linearen Funktion entsprechen, für die ich zwei lineare Regressionen ausführen möchte. Es gibt einen Punkt, an dem sich die Steigung und der Achsenabschnitt ändern, und ich muss diesen Punkt finden (ein Programm schreiben, um ihn zu finden).

Gedanken?

Das mcpPaket kann dies tun. Angenommen, Ihre Daten sind

Lassen Sie uns zunächst einige Daten simulieren:

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

Nun wollen wir sehen, ob wir den Änderungspunkt bei 40 (und die Parameterwerte) wiederherstellen können, indem wir mcp:

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

Plotten Sie es. Die grauen Linien sind zufällige Zeichnungen aus der Anpassung, die zeigen, dass sie den Trend erfassen. Die blaue Kurve ist die geschätzte Position des Änderungspunkts:

Sehen wir uns die einzelnen Parameterschätzungen an. int_sind Abschnitte, x_sind Steigungen auf x und cp_sind Änderungspunkte:

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

Haftungsausschluss: Ich bin der Entwickler von mcp.

R-Paket strucchange könnte Ihnen helfen. Schauen Sie sich die Vignette an, sie hat einen schönen Überblick darüber, wie man ähnliche Probleme löst.

$X_i=(x_i,y_i)$$i=1,..,N$$j$$2$$N-2$$\{X_1,...,X_j\}$$\{X_{(j+1)},...,X_N\}$$j$

Dies ist ein (Offline-) Änderungspunkterkennungsproblem. Unsere vorherige Diskussion enthält Verweise auf Zeitschriftenartikel und R-Code. Schauen Sie sich zuerst das "Produktpartitionsmodell" von Barry und Hartigan an , da es Änderungen in der Steigung handhabt und effizient implementiert werden kann.

Auch das segmentierte Paket hat mir in der Vergangenheit bei ähnlichen Problemen geholfen.

Ich baute auf der Antwort von mbq auf, nach allen Möglichkeiten zu suchen. Außerdem mache ich das:

• Überprüfen Sie die Signifikanz der beiden stückweisen Modelle, um sicherzustellen, dass die Koeffizienten signifikant sind
• Überprüfen Sie die Differenz zur Summe der quadratischen Residuen für das vollständige Modell
• Bestätige mein Modell visuell (stelle sicher, dass es kein Unsinn ist)

Warum auf die Bedeutung prüfen? Dies liegt daran, dass der Punkt mit der minimalen SSE bedeutungslos ist, wenn eines der stückweisen Modelle sehr schlecht zu den Daten passt. Dies kann für zwei stark korrelierte Variablen ohne einen klaren Haltepunkt geschehen, an dem sich die Steigungen ändern.

Lassen Sie uns diesen einfachen Ansatz anhand eines einfachen Testfalls überprüfen:

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

Der Haltepunkt ist offensichtlich Null. Verwenden Sie das folgende R-Skript:

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
}

Passen Sie stückweise lineare Modelle für alle möglichen Kombinationen an:

f(x,y)
   k sums
   0    0

Wenn wir die Koeffizienten für die beiden optimalen Modelle überprüfen, sind sie von hoher Bedeutung. Ihr R2 wird auch sehr hoch sein.