R: Интерполировать + Экстраполировать через границы - с data.table?

Итак, у меня есть следующая проблема: у меня есть набор данных A (объект data.table) следующей структуры:

date days rate 1996-01-02 9 5.763067 1996-01-02 15 5.745902 1996-01-02 50 5.673317 1996-01-02 78 5.608884 1996-01-02 169 5.473762 1996-01-03 9 5.763067 1996-01-03 14 5.747397 1996-01-03 49 5.672263 1996-01-03 77 5.603705 1996-01-03 168 5.470584 1996-01-04 11 5.729460 1996-01-04 13 5.726104 1996-01-04 48 5.664931 1996-01-04 76 5.601891 1996-01-04 167 5.468961

Обратите внимание, что столбец days и его размер могут различаться для каждого дня. Моя цель теперь состоит в том, чтобы (кусочно-линейно) интерполировать скорость по дням. Я делаю это каждый день через

approx(x=A[,days],y=A[,rate],xout=days_vec,rule=2)

где days_vec <- min_days:max_days, то есть интересующий меня диапазон дней (скажем, 1:100).

У меня тут две проблемы:

  1. приблизительно только интерполирует, т. е. не создает линейную подгонку по min(x) и max(x). Если меня сейчас интересуют дни 1:100, мне сначала нужно сделать это вручную, используя дни 9 и 15 (первые 2 строки A) через:

    first_days <- 1:(A[1,days]-1) #1:8 rate_vec[first_days] <- A[1,rate] + (first_days - A[1,days])/(A[2,days]-A[1,days])*(A[2,rate]-A[1,rate])

а затем используйте строку примерно выше для rate_vec[9:100]. Есть ли способ сделать это за 1 шаг?

  1. Прямо сейчас, учитывая, что мне нужно два шага, а точка сдвига между двумя процедурами (здесь 9) различается в зависимости от даты, я не вижу реализации через data.table, хотя это было бы гораздо предпочтительнее (используя методы data.table для интерполяции/ экстраполировать, а затем вернуть расширенный объект data.table). Таким образом, в настоящее время я запускаю цикл for по датам, что, конечно, намного медленнее.

Вопрос. Можно ли решить описанную выше проблему лучше, а также можно ли это каким-то образом решить с помощью методов data.table вместо цикла A?


r data.table interpolation extrapolation
person Daedalus    schedule 05.09.2017    source источник
comment
Можете ли вы быть немного более ясным о том, что вы пытаетесь сделать? Похоже, вы пытаетесь сделать три вещи: (1) создать дополнительные строки для значений days, которых нет в таблице (2) создать lm для rate ~ days для каждой группы (определяется date) (3) использовать predict.lm для создания значения rate для days, где rate отсутствует   -  person C8H10N4O2    schedule 05.09.2017

Ответы (1)


arrow_upward
5
arrow_downward

Как насчет чего-то подобного.

# please try to make a fully reproducible example!
library(data.table)
df <- fread(input=
"date       days     rate
1996-01-02    9 5.763067
1996-01-02   15 5.745902
1996-01-02   50 5.673317
1996-01-02   78 5.608884
1996-01-02  169 5.473762
1996-01-03    9 5.763067
1996-01-03   14 5.747397
1996-01-03   49 5.672263
1996-01-03   77 5.603705
1996-01-03  168 5.470584
1996-01-04   11 5.729460
1996-01-04   13 5.726104
1996-01-04   48 5.664931
1996-01-04   76 5.601891
1996-01-04  167 5.468961")
df[,date := as.Date(date)]

1. Создайте значения скорости NA для дней в диапазоне 1: 100, которых нет в наборе данных.

df <- 
  merge(df,
        expand.grid( days=1L:100L,                   # whatever range you are interested in
                     date=df[,sort(unique(date))] ), # dates with at least one observation
        all=TRUE # "outer join" on all common columns (date, days)
        )

2. Для каждого значения даты используйте линейную модель, чтобы предсказать значения скорости NA.

df[, rate := ifelse(is.na(rate),
                    predict(lm(rate~days,.SD),.SD), # impute NA w/ lm using available data
                    rate),                          # if not NA, don't impute
   keyby=date]

Дает тебе:

head(df,10)
#           date days     rate
#  1: 1996-01-02    1 5.766787 <- rates for days 1-8 & 10 are imputed 
#  2: 1996-01-02    2 5.764987
#  3: 1996-01-02    3 5.763186
#  4: 1996-01-02    4 5.761385
#  5: 1996-01-02    5 5.759585
#  6: 1996-01-02    6 5.757784
#  7: 1996-01-02    7 5.755983
#  8: 1996-01-02    8 5.754183
#  9: 1996-01-02    9 5.763067 <- this rate was given
# 10: 1996-01-02   10 5.750581

Если есть значения date без хотя бы двух наблюдений rate, вы, вероятно, получите сообщение об ошибке, потому что у вас не будет достаточно точек для размещения линии.

Альтернатива: решение для прокатных соединений для кусочно-линейной интерполяции

Для этого требуется скользящее соединение слева и справа, а также среднее значение двух, которое игнорирует значения NA.

Однако это не очень хорошо для экстраполяции, поскольку это просто константа (либо первое, либо последнее наблюдение) вне индексов наблюдений.

setkey(df, date, days)

df2 <- data.table( # this is your framework of date/days pairs you want to evaluate
        expand.grid( date=df[,sort(unique(date))],
                     days=1L:100L),
        key = c('date','days')
)

# average of non-NA values between two vectors
meanIfNotNA <- function(x,y){
  (ifelse(is.na(x),0,x) + ifelse(is.na(y),0,y)) /
    ( as.numeric(!is.na(x)) + as.numeric(!is.na(y)))
}

df3 <- # this is your evaluations for the date/days pairs in df2.
  setnames(
    df[setnames( df[df2, roll=+Inf], # rolling join Last Obs Carried Fwd (LOCF)
                 old = 'rate',
                 new = 'rate_locf' 
               ),
     roll=-Inf], # rolling join Next Obs Carried Backwd (NOCB)
    old = 'rate',
    new = 'rate_nocb'
  )[, rate := meanIfNotNA(rate_locf,rate_nocb)]
# once you're satisfied that this works, you can include rate_locf := NULL, etc.

head(df3,10)
#          date days rate_nocb rate_locf     rate
#  1: 1996-01-02    1  5.763067        NA 5.763067
#  2: 1996-01-02    2  5.763067        NA 5.763067
#  3: 1996-01-02    3  5.763067        NA 5.763067
#  4: 1996-01-02    4  5.763067        NA 5.763067
#  5: 1996-01-02    5  5.763067        NA 5.763067
#  6: 1996-01-02    6  5.763067        NA 5.763067
#  7: 1996-01-02    7  5.763067        NA 5.763067
#  8: 1996-01-02    8  5.763067        NA 5.763067
#  9: 1996-01-02    9  5.763067  5.763067 5.763067  <- this rate was given
# 10: 1996-01-02   10  5.745902  5.763067 5.754485
person C8H10N4O2    schedule 05.09.2017
comment
Я не видел кусочно-линейного требования выше - если это не работает для вас, и вы действительно просто хотите получить информацию из следующих наблюдений слева и справа, я могу отредактировать, чтобы включить скользящие соединения. Но сначала объясните, как вы хотите обрабатывать конечные точки (т. е. пропущенные значения до первого наблюдения и после последнего наблюдения). - person C8H10N4O2; 05.09.2017
comment
Большое спасибо!! Это определенно очень помогает и работает очень гладко. Причина, по которой я хотел, чтобы она была кусочно-линейной, заключается в том, что скорости следуют не линейной схеме, а скорее вогнутой, которую я хотел аппроксимировать несколькими кусочно-линейными интерполяциями, чтобы приблизиться к истинной функции. Моя идея заключалась в том, чтобы обрабатывать пропущенные значения, продолжая линейный тренд между двумя предыдущими узлами (см. функцию линейной интерполяции для случая 1:8 выше). К сожалению, я еще не очень хорошо знаком со скользящими соединениями пакета data.table, поэтому, если бы вы могли отредактировать это выше, было бы здорово! - person Daedalus; 05.09.2017
comment
Альтернативой может быть использование всей кривой с вашим подходом и использование сплайнов вместо линейных моделей и прогнозирование значений na на основе сглаженных сплайнов. Я также рассмотрю это. - person Daedalus; 05.09.2017
comment
@ Daedalus, вот подвижные соединения. Обратите внимание, что он использует две ближайшие точки для линейной интерполяции и только одну для экстраполяции, которая просто дает вам константу. В зависимости от того, насколько вы заботитесь об экстраполяции, вам может быть лучше использовать сплайны или какой-то действительно индивидуальный подход. Вы также можете комбинировать интерполяцию с скользящими соединениями, а затем что-то еще для экстраполяции. Цель этого ответа состояла в том, чтобы помочь вам эффективно работать в data.table. - person C8H10N4O2; 05.09.2017
comment
Спасибо, кофеинщик! Это очень помогло. Я думаю, что попробую как сплайн-версию с data.table (которая должна быть реализована так же, как и линейные модели, как я думаю), так и линейные прокатные соединения, и сравним, что лучше подходит. Ваше здоровье - person Daedalus; 05.09.2017