r - 评估公式的快速方法？

Question

我使用 dyn 或 dynlm 来预测使用滞后变量的时间序列。

然而，在这两种情况下，预测函数一次只计算一个时间步，在我的计算机上每步需要 24 毫秒的恒定时间，或者我的数据集需要大约 1.8 小时，这是超长的，因为整个回归需要大约10 秒。

所以，我在想，也许最快的方法可能只是手动评估公式？

那么，是否有某种方法可以评估给定 data.frame 或当前环境或类似环境中给定值的公式？

我正在考虑以下方面的事情:

evalMagic( load ~ temperature + time, data.frame( temperature = 10, time = 4 ) )

我想，在我写这篇文章的时候，我们需要以某种方式处理这些系数，比如:

evalMagic( load ~ temperature + time, data.frame( temperature = 10, time = 4 ), model$coefficients )

.... 所以这就提出了以下问题:

这不是 predict 应该做的吗？

为什么预测这么慢？

我有什么选择可以使预测更快一点？毕竟，它不是反转任何矩阵或其他东西，它只是一点算术!

Answer 1

最后我写了自己的滞后实现。它很笨拙，也不漂亮，但速度要快得多。它可以在我蹩脚的笔记本电脑上在 4 秒内处理 1000 行。

# lags is a data.frame, eg:
#   var  amount
#   y    1
#   y    2
addLags <- function( dataset, lags ) {
    N <- nrow(dataset)
    print(lags)
    if( nrow(lags) > 0 ) {
        print(lags)
        for( j in 1:nrow(lags) ) {
            sourcename <- as.character( lags[j,"var"] )
            k <- lags[j,"amount"]
            cat("k",k,"sourcename",sourcename,"\n")
            lagcolname <- sprintf("%s_%d",sourcename,k)
            dataset[,lagcolname] <- c(rep(0,k), dataset[1:(N-k),sourcename])
        }
    }
    dataset
}

lmLagged <- function( formula, train, lags ) {
    # get largest lag, and skip that
    N <- nrow(train)
    skip <- 0
    for( j in 1:nrow(lags) ) {
        k <- lags[j,"amount"]
        skip <- max(k,skip)
    }
    print(train)
    train <- addLags( train, lags )
    print(train)
    lm( formula, train[(skip+1):N,] )
}

# pass in training data, test data,
# it will step through one by one
# need to give dependent var name
# lags is a data.frame, eg:
#   var amount
#   y    1
#   y    2
predictLagged <- function( model, train, test, dependentvarname, lags ) {
    Ntrain <- nrow(train)
    Ntest <- nrow(test)
    test[,dependentvarname] <- NA
    testtraindata <- rbind( train, test )
    testtraindata <- addLags( testtraindata, lags )
    for( i in 1:Ntest ) {
       thistestdata <- testtraindata[Ntrain + i,]
       result <- predict(model,newdata=thistestdata)
       for( j in 1:nrow(lags) ) {
            sourcename <- lags[j,"var"]
            k <- lags[j,"amount"]
            lagcolname <- sprintf("%s_%d",sourcename,k)
            testtraindata[Ntrain + i + k,lagcolname] <- result
       }
       testtraindata[Ntrain+i,dependentvarname] <- result
    }
    return( testtraindata[(Ntrain+1):(Ntrain + Ntest),dependentvarname] )    
}

library("RUnit")

# size of training data
N <- 6
predictN <- 50

# create training data, which we can get exact fit on
set.seed(1)
x = sample( 100, N )
traindata <- numeric()
traindata[1] <- 1 + 1.1 * x[1]
traindata[2] <- 2 + 1.1 * x[2]
for( i in 3:N ) {
   traindata[i] <- 0.5 + 0.3 * traindata[i-2] - 0.8 * traindata[i-1] + 1.1 * x[i]
}
train <- data.frame(x = x, y = traindata, foo = 1)
#train$x <- NULL

# create testing data, bunch of NAs
test <- data.frame( x = sample(100,predictN), y = rep(NA,predictN), foo = 1)

# specify which lags we need to handle
# one row per lag, with name of variable we are lagging, and the distance
# we can then use these in the formula, eg y_1, and y_2
# are y lagged by 1 and 2 respectively
# It's hacky but it kind of works...
lags <- data.frame( var = c("y","y"), amount = c(1,2) ) 

# fit a model
model <- lmLagged(  y ~ x + y_1 + y_2, train, lags )
# look at the model, it's a perfect fit. Nice!
print(model)

print(system.time( test <- predictLagged( model, train, test, "y", lags ) ))
#checkEqualsNumeric( 69.10228, test[56-6], tolerance = 0.0001 )
#checkEquals( 2972.159, test$y[106-6] )
print(test)

# nice plot
plot(test, type='l')

输出:

> source("test/test.regressionlagged.r",echo=F)

Call:
lm(formula = formula, data = train[(skip + 1):N, ])

Coefficients:
(Intercept)            x          y_1          y_2  
        0.5          1.1         -0.8          0.3  

   user  system elapsed 
  0.204   0.000   0.204 
 [1]  -19.108620  131.494916  -42.228519   80.331290  -54.433588   86.846257
 [7]  -13.807082   77.199543   12.698241   64.101270   56.428457   72.487616
[13]   -3.161555   99.575529    8.991110   44.079771   28.433517    3.077118
[19]   30.768361   12.008447    2.323751   36.343533   67.822299  -13.154779
[25]   72.070513  -11.602844  115.003429  -79.583596  164.667906 -102.309403
[31]  193.347894 -176.071136  254.361277 -225.010363  349.216673 -299.076448
[37]  400.626160 -371.223862  453.966938 -420.140709  560.802649 -542.284332
[43]  701.568260 -679.439907  839.222404 -773.509895  897.474637 -935.232679
[49] 1022.328534 -991.232631

这 91 行代码大约需要 12 个小时的工作时间。好吧，我承认我玩了一段时间植物和僵尸。所以，10个小时。加上午餐和晚餐。尽管如此，还是有很多工作要做。

如果我们将 predictN 更改为 1000，我会从 system.time 调用中得到大约 4.1 秒。

我认为它更快，因为:

我们不使用时间序列；我怀疑这会加快速度

我们不使用动态lm库，只使用普通lm；我想这会稍微快一点

我们只将单行数据传递给每个预测的预测，我认为这要快得多，例如使用 dyn$lm 或 dynmlm，如果滞后为 30，则需要将 31 行数据传递给预测 AFAIK

少了很多 data.frame/matrix 复制，因为我们只是在每次迭代时就地更新滞后值

编辑:更正了小错误，其中 predictLagged 返回了一个多列数据帧，而不仅仅是一个数字向量
编辑 2:更正了不能添加多个变量的较小错误。还协调了滞后的注释和代码，并将滞后结构更改为“var”和“amount”以代替“name”和“lags”。此外，更新了测试代码以添加第二个变量。

编辑:这个版本有很多错误，我知道，因为我对它进行了更多的单元测试并修复了它们，但是复制和粘贴非常耗时，所以我会在几天后更新这篇文章，一旦我的截止日期结束。