Оптимизированные функции качения на нерегулярных временных рядах с временным окном
Есть ли способ использовать rollapply (из zoo
package или что-то подобное) оптимизированные функции (rollmean
, rollmedian
и т.д.) для вычисления функций качения с использованием временного окна вместо одного, основанного на числе наблюдений? Я хочу просто: для каждого элемента в нерегулярном временном ряду я хочу вычислить функцию качения с окном N дней. То есть, окно должно включать все наблюдения за N дней до текущего наблюдения. Временные ряды также могут содержать дубликаты.
Здесь следует пример. Учитывая следующие временные ряды:
date value
1/11/2011 5
1/11/2011 4
1/11/2011 2
8/11/2011 1
13/11/2011 0
14/11/2011 0
15/11/2011 0
18/11/2011 1
21/11/2011 4
5/12/2011 3
Скользящая медиана с 5-дневным окном, выровненным вправо, должна приводить к следующему вычислению:
> c(
median(c(5)),
median(c(5,4)),
median(c(5,4,2)),
median(c(1)),
median(c(1,0)),
median(c(0,0)),
median(c(0,0,0)),
median(c(0,0,0,1)),
median(c(1,4)),
median(c(3))
)
[1] 5.0 4.5 4.0 1.0 0.5 0.0 0.0 0.0 2.5 3.0
Я уже нашел некоторые решения там, но они обычно сложны, что обычно означает медленное. Мне удалось реализовать собственный калькулятор функции. Проблема в том, что для очень длинных рядов оптимизированная версия медианного (rollmedian) может сделать огромную разницу во времени, поскольку она учитывает перекрытие между окнами. Я бы хотел избежать повторного его реализации. Я подозреваю, что есть некоторые трюки с параметрами rollapply, которые заставят его работать, но я не могу понять это. Заранее спасибо за помощь.
Ответы
Ответ 1
В большинстве ответов предлагается вставить NA, чтобы временные ряды были регулярными.
Однако это может быть медленным в случае длинных временных рядов. Кроме того, он не работает для функций, которые нельзя использовать с NA.
Аргумент ширины rollapply (пакет zoo) может быть списком (подробнее см. справку rollapply). На основании этого я написал функцию, которая создает список, который будет использоваться с rollapply в качестве параметра ширины. Функция извлекает индексы для нерегулярных объектов зоопарка, если движущееся окно должно быть временным, а не индексированным. Поэтому индекс объекта зоопарка должен быть фактическим.
# Create a zoo object where index represents time (e.g. in seconds)
d <- zoo(c(1,1,1,1,1,2,2,2,2,2,16,25,27,27,27,27,27,31),
c(1:5,11:15,16,25:30,31))
# Create function
createRollapplyWidth = function(zoodata, steps, window ){
mintime = min(time(zoodata))
maxtime = max(time(zoodata))
spotstime = seq(from = mintime , to = maxtime, by = steps)
spotsindex = list()
for (i in 1:length(spotstime)){
spotsindex[[i]] = as.numeric(which(spotstime[i] <= time(zoodata) & time(zoodata) < spotstime[i] + window))}
rollapplywidth = list()
for (i in 1:length(spotsindex)){
if (!is.na(median(spotsindex[[i]])) ){
rollapplywidth[[round(median(spotsindex[[i]]))]] = spotsindex[[i]] - round(median(spotsindex[[i]]))}
}
return(rollapplywidth)
}
# Create width parameter for rollapply using function
rollwidth = createRollapplyWidth(zoodata = d, steps = 5, window = 5)
# Use parameter in rollapply
result = rollapply(d, width = rollwidth , FUN = sum, na.rm = T)
result
Ограничение: не основано на дате, но в секундах. Параметр "частичный" rollapply не работает.
Ответ 2
Вот моя работа с проблемой. Если такой подход зависит от того, что вы хотели (я не знаю, удовлетворительно ли это с точки зрения скорости), я могу написать его как более подробный ответ (хотя он основан на идее @rbatt).
library(zoo)
library(dplyr)
# create a long time series
start <- as.Date("1800-01-01")
end <- as.Date(Sys.Date())
df <- data.frame(V1 = seq.Date(start, end, by = "day"))
df$V2 <- sample(1:10, nrow(df), replace = T)
# make it an irregular time series by sampling 10000 rows
# including allowing for duplicates (replace = T)
df2 <- df %>%
sample_n(10000, replace = T)
# create 'complete' time series & join the data & compute the rolling median
df_rollmed <- data.frame(V1 = seq.Date(min(df$V1), max(df$V1), by = "day")) %>%
left_join(., df2) %>%
mutate(rollmed = rollapply(V2, 5, median, na.rm = T, align = "right", partial = T)) %>%
filter(!is.na(V2)) # throw out the NAs from the complete dataset
Ответ 3
Не проверяйте скорость, но если дата не имеет более чем max.dup
, то должно быть, что последние 5 * max.dup записей содержат последние 5 дней, поэтому приведенная ниже однострочная функция fn
на rollapplyr
сделает это:
k <- 5
dates <- as.numeric(DF$date)
values <- DF$value
max.dup <- max(table(dates))
fn <- function(ix, d = dates[ix], v = values[ix], n = length(ix)) median(v[d >= d[n]-k])
rollapplyr(1:nrow(DF), max.dup * k, fn, partial = TRUE)
## [1] 5.0 4.5 4.0 1.0 0.5 0.0 0.0 0.0 2.5 3.0
Примечание: Мы использовали это для DF
:
Lines <- "
date value
1/11/2011 5
1/11/2011 4
1/11/2011 2
8/11/2011 1
13/11/2011 0
14/11/2011 0
15/11/2011 0
18/11/2011 1
21/11/2011 4
5/12/2011 3
"
DF <- read.table(text = Lines, header = TRUE)
DF$date <- as.Date(DF$date, format = "%d/%m/%Y")
Ответ 4
Мы можем сделать это, просто используя базу, следующим образом:
Сначала настройте данные (на основе примечания @g-grothendieck)
library(data.table)
Lines <- "
date value
1/11/2011 5
1/11/2011 4
1/11/2011 2
8/11/2011 1
13/11/2011 0
14/11/2011 0
15/11/2011 0
18/11/2011 1
21/11/2011 4
5/12/2011 3
"
DT <- as.data.table(read.table(text = Lines, header = TRUE))
DT$date <- as.Date(DF$date, format = "%d/%m/%Y")
DT$row <- 1:NROW(DF)
setkey(DT, row, date) #mark columns as sorted, for speed
Обратите внимание, что я добавил вектор в таблицу данных, содержащую номер строки, чтобы мы могли передать номер строки в функцию apply. Я также использовал таблицу данных, чтобы упростить синтаксис для следующего шага и ускорить эту функцию, если она применяется к большим массивам. Теперь мы применяем следующее:
roll.median.DT <- function(x){
this.date <- as.Date(x[1])
this.row <- as.numeric(x[3])
median(DT[row <= this.row & date >= (this.date-5)]$value) #NB DT is not defined within function, so it is found from parent scope
}
apply(DT, FUN=roll.median.DT, MARGIN = 1)
[1] 5.0 4.5 4.0 1.0 0.5 0.0 0.0 0.0 2.5 3.0