Оптимизация кусочно-постоянной подгонки в R для многих фреймов данных - PullRequest
Новая
       72

Оптимизация кусочно-постоянной подгонки в R для многих фреймов данных

0 голосов
/ 29 октября 2018

Я хочу оптимизировать код, чтобы соответствовать постоянной кусочной функции для моих данных. У меня есть df, как это: 

df = data.frame (x = 1:180,
                 y = c(0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,1,0,
                       0,0,0,0,0,0,0,0,0,0,1,2,0,0,0,2,2,4,2,2,3,2,1,2,0,1,0,1,4,
                       0,1,2,3,1,1,1,0,2,0,3,2,1,1,1,1,5,4,2,1,0,2,1,1,2,0,0,2,2,
                       1,1,1,0,0,0,0,2,3,0,3,2,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
                       0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
                       0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,0,
                       0,0,0,0,0,0))

Со вторым ответом здесь , который автор назвал "грубой силой", я написал следующий код, чтобы применить тот же метод к 960 фреймам данных: 

# Set the possible values of k1 and k2
k1 <- df$x[df$x < df$x[df$y == max(df$y)]]
k2 <-  df$x[df$x > df$x[df$y == max(df$y)]]
# get all combinations of k1 and k2
grid <- subset(expand.grid(k1 = k1, k2 = k2), k1 < k2)

store <- numeric(nrow(grid))

for(j in 1:nrow(grid)){
  k1 <- grid$k1[j]
  k2 <- grid$k2[j]
  model <- lm( y ~ I(x < k1) + I(x >= k1 & x < k2) + I(x >= k2),data = df)
# store sigmas for every pair of values
  store[j] <-  sigma(model)
}
# Get those values which minimize the residual standard deviation
bp <- grid[which.min(store),]

# Look for the minimum sigma (dot in red)
plot(unlist(store),pch =21,bg = "lightblue",
     ylab = "Sigma",xlab = "Pairs of k1 and k2")
points(which.min(store),
       unlist(store)[which.min(store)],col = "red",pch = 16)

enter image description here

# Fit a model with the best k

fit <- lm(y ~ I(x < bp[1,"k1"]) + I(x >= bp[1,"k1"] & x < bp[1,"k2"]) + I(x >= bp[1,"k2"]),data = df)

plot(df,pch = 21,bg = "orange")
lines(fitted(fit),col = "blue",lwd = 2)

enter image description here

Этот код работает нормально (я провожу некоторые визуальные проверки в случайных кадрах данных), но на моем компьютере это занимает почти 6 часов. Может ли кто-нибудь предложить лучший способ достичь того же результата? Если это можно сделать, поддерживая читаемый код, это будет лучше, но это не обязательно.

Заранее спасибо

...