R：在边界上插值 + 外插 - 使用 data.table？
R: Interpolate + Extrapolate over borders - with 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
请注意，天 列及其大小可能每天都不同。
我现在的目标是（分段线性）沿天插值率。我每天 通过
这样做
approx(x=A[,days],y=A[,rate],xout=days_vec,rule=2)
其中 days_vec <- min_days:max_days，即我感兴趣的天数范围（比如 1:100）。
我这里有两个问题：
approx 仅进行插值，即它不会在 min(x) 和 max(x) 之间创建线性拟合。如果我现在对 1:100 天感兴趣，我首先需要使用第 9 天和第 15 天（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])
然后使用上面的 approx 行 rate_vec[9:100]。有没有一种方法可以一步完成？
- 现在，考虑到我需要两个步骤并且两个过程之间的转换点（此处为 9）因日期而异，我看不到通过 data.table 的实现，尽管这是非常可取的（使用data.table 方法到 interpolate/extrapolate，然后返回扩展的 data.table 对象）。因此，我目前 运行 一个 for 循环遍历日期，这当然要慢得多。
问题：上面的问题是否更容易实现，而且，这是否可以通过 data.table 方法而不是通过 A 循环来实现？
这样的怎么样。
# 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。为 1:100 范围内但不在数据集
中的天数创建 NA 值
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
所以，我遇到了以下问题： 我有一个数据集 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
请注意，天 列及其大小可能每天都不同。 我现在的目标是（分段线性）沿天插值率。我每天 通过
这样做approx(x=A[,days],y=A[,rate],xout=days_vec,rule=2)
其中 days_vec <- min_days:max_days，即我感兴趣的天数范围（比如 1:100）。
我这里有两个问题：
approx 仅进行插值，即它不会在 min(x) 和 max(x) 之间创建线性拟合。如果我现在对 1:100 天感兴趣，我首先需要使用第 9 天和第 15 天（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])
然后使用上面的 approx 行 rate_vec[9:100]。有没有一种方法可以一步完成？
- 现在，考虑到我需要两个步骤并且两个过程之间的转换点（此处为 9）因日期而异，我看不到通过 data.table 的实现，尽管这是非常可取的（使用data.table 方法到 interpolate/extrapolate，然后返回扩展的 data.table 对象）。因此，我目前 运行 一个 for 循环遍历日期，这当然要慢得多。
问题：上面的问题是否更容易实现，而且，这是否可以通过 data.table 方法而不是通过 A 循环来实现？
这样的怎么样。
# 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。为 1:100 范围内但不在数据集
中的天数创建 NA 值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