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

10

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?

Rhombidodekaeder
quelle
3
Was ist die Richtlinie für Cross-Posting? Die exakt gleiche Frage wurde auf math.stackexchange.com gestellt: math.stackexchange.com/questions/15214/…
mpiktas
Was ist in diesem Fall falsch daran, einfache nichtlineare kleinste Quadrate zu machen? Vermisse ich etwas Offensichtliches?
Grg s
Ich würde sagen, dass die Ableitung der Zielfunktion in Bezug auf den Änderungspunktparameter ziemlich ungleichmäßig ist
Andre Holzner
Die Steigung würde sich so stark ändern, dass ein nichtlineares kleinstes Quadrat nicht präzise und genau wäre. Was wir wissen ist, dass wir zwei oder mehr lineare Modelle haben, daher sollten wir zuschlagen, um diese beiden Modelle zu extrahieren.
HelloWorld

Antworten:

1

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:

Geben Sie hier die Bildbeschreibung ein

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.

Jonas Lindeløv
quelle
8

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.

mpiktas
quelle
6

X.ich=(xich,yich)ich=1,..,N.j2N.- -2{X.1,...,X.j}}{X.(j+1),...,X.N.}}j


quelle
Ich habe eine Antwort auf der Grundlage Ihres einfachen, aber effektiven Vorschlags veröffentlicht.
HelloWorld
5

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.

whuber
quelle
3

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

Mischa
quelle
Leider benötigt das Paket einen Startwert für den Haltepunkt.
HelloWorld
Außerdem segmentedkönnen Intercept-Änderungen zwischen Segmenten nicht modelliert werden - nur ein Intercept für das erste Segment.
Jonas Lindeløv
2

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)

Geben Sie hier die Bildbeschreibung ein

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.

Hallo Welt
quelle