Menemukan titik perubahan dalam data dari fungsi linear piecewise

Salam pembuka,

Saya melakukan penelitian yang akan membantu menentukan ukuran ruang yang diamati dan waktu yang berlalu sejak big bang. Semoga Anda bisa membantu!

Saya memiliki data yang sesuai dengan fungsi linear piecewise di mana saya ingin melakukan dua regresi linier. Ada titik di mana perubahan kemiringan dan mencegat, dan saya perlu (menulis program untuk) menemukan titik ini.

Pikiran?

The mcppaket dapat melakukan hal ini. Katakan data Anda

Pertama, mari kita simulasikan beberapa data:

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

Sekarang mari kita lihat apakah kita dapat memulihkan titik perubahan di 40 (dan nilai parameter) menggunakan mcp:

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

Plot itu. Garis abu-abu adalah undian acak dari fit, menunjukkan bahwa itu menangkap tren. Kurva biru adalah perkiraan lokasi titik perubahan:

Mari kita lihat estimasi parameter individual. int_adalah penyadapan, x_adalah kemiringan pada x, dan cp_merupakan titik perubahan:

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

Penafian: Saya adalah pengembang dari mcp.

Strucchange paket R dapat membantu Anda. Lihatlah sketsa, ia memiliki ikhtisar yang bagus bagaimana menyelesaikan masalah yang sama.

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

Ini adalah masalah pendeteksian changepoint (offline). Diskusi kami sebelumnya memberikan referensi ke artikel jurnal dan kode R. Lihatlah dulu "model partisi produk" Barry dan Hartigan , karena menangani perubahan kemiringan dan memiliki implementasi yang efisien.

Juga tersegmentasi paket telah membantu saya dengan masalah serupa di masa lalu.

Saya membangun berdasarkan jawaban mbq yang mencari semua kemungkinan. Selanjutnya, saya melakukan ini:

• Periksa signifikansi dari kedua model sambungan untuk memastikan koefisiennya signifikan
• Periksa perbedaan dengan jumlah residu kuadrat untuk model lengkap
• Konfirmasikan model saya secara visual (pastikan itu bukan sesuatu yang tidak masuk akal)

Mengapa memeriksa signifikansi? Itu karena titik dengan SSE minimum tidak ada artinya jika salah satu model sambungan sesuai data yang sangat buruk. Hal ini dapat terjadi untuk dua variabel yang sangat berkorelasi tanpa breakpoint yang jelas di mana perubahan lereng.

Mari kita periksa pendekatan sederhana ini dengan test case yang mudah:

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

Breakpoint jelas nol. Gunakan skrip R berikut:

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
}

Sesuaikan model linear satu demi satu untuk semua kemungkinan kombinasi:

f(x,y)
   k sums
   0    0

Jika kita periksa koefisien untuk dua model optimal, mereka akan sangat signifikan. R2 mereka juga akan sangat tinggi.