創建遞歸式樣本外估計以有效計算R中的RMSE的方法

我試圖評估不同OLS模型的樣本外預測性能。最簡單的時間序列迴歸看起來是這樣的： Y_t = B0 + B1 * Y_t-30 + E_T創建遞歸式樣本外估計以有效計算R中的RMSE的方法

配件期爲模型，比方說50，然後我讓模型運行使用dynlm包

dynlm(as.zoo(Y) ~ L(as.zoo(Y), 30), start = "1996-01-01", end = timelist[i])

在我當前的代碼中，我只是讓索引運行到結束，然後保存相應模型的RMSE。但是這個RMSE並不是預測中的一步，因爲我現在的代碼已經很慢了，甚至沒有做到我想做的事，所以我想問你是否有一個建議我應該用它來實現我的目標。

概括起來講，我要做到以下幾點：

1）一定擬合時間後運行一個遞推回歸（擴大窗口，而不是滾動窗口）

2）提前創建的一步出的樣本外預測

3）計算均方根這些預測與實際觀測的方差來評估模型的性能

我試着這樣做，到目前爲止有一個巨大的for循環和dynlm包，但水庫ults並不是真的令人滿意。任何輸入是非常感激，因爲我一直在尋找解決方案相當長一段時間了。我會盡快更新我的示例代碼。

# minimal working example 
require(xts) 
require(zoo) 
require(dynlm) 
timelist <- seq.Date(from = as.Date("1996-01-01"), to = as.Date("1998-01-01"), by = "days") 
set.seed(123) 
Y <- xts(rnorm(n = length(timelist)), order.by = timelist) 
X <- xts(rnorm(n = length(timelist), mean = 10), order.by = timelist) 
# rmse container 
rmse.container.full <- data.frame(matrix(NA, ncol = 3, nrow = length(index(timelist)))) 
colnames(rmse.container.full) <- c("Date", "i", "rmse.m1") 
rmse.container.full$Date <- timelist 
# fitting period 
for(i in 50:length(timelist)) { 
    # m1 
    model1 <- dynlm(as.zoo(Y) ~ L(as.zoo(X), 30), start = "1996-01-01", end = timelist[i]) 
    rmse.container.full[i, 2] <- i 
    rmse.container.full[i, 3] <- summary(model1)$sigma # RSME mod1 etc 
    print(i) 
}

來源

2016-06-24 tester

歡迎StackOverflow上這樣做。請看看這些關於如何產生[最小，完整和可驗證的例子]（http://stackoverflow.com/help/mcve）的技巧，以及這篇文章[在R中創建一個很好的例子]（ http://stackoverflow.com/questions/5963269/how-to-make-a-great-r-reproducible-example）。也許下面的提示[問一個好問題]（http://stackoverflow.com/help/how-to-ask）也值得一讀。 – lmo

爲什麼不是結果「非常令人滿意」？這不是我們可以解決的具體問題。如果你需要統計建模方面的幫助，你應該在[stats.se]中提問你的問題，否則請清楚這裏的編程問題是什麼，並且包括一個[可重現的例子]（http://stackoverflow.com/questions/5963269/ how-to-make-a-great-r-reproducible-example） – MrFlick

感謝您的建議@MrFlick，我現在看看這個： [鏈接]（https://www.otexts.org/fpp ）並檢查交叉驗證。 – tester

那麼，作爲一個誰問我願意貢獻我如何解決我的問題，這樣的問題：

正如我只需要領先一步的預測，我可以丟掉一切，本作的代碼運行速度更快。（從12分鐘到每個型號約10秒）。

我自己創建了完整的數據框（包括滯後），並使用lm而不是dynlm。下面的代碼給了我我想要的結果（我手動檢查了前幾個觀察結果，它似乎工作）。該代碼是改編自點擊這裏：Recursive regression in R

 mod1.predictions <- lapply(seq(1400, nrow(df.full)-1), function(x) { 
       mod1 <- lm(Y ~ X, data = df.full[1:x, ]) 
       pred1 <- predict(mod1, newdata = df.full[x+1, ]) 
       return(pred1) 
       })

爲了計算RMSE我用這個功能

# rmse function 
rmse <- function(sim, obs) { 
    res <- sqrt(mean((sim - obs)^2, na.rm = TRUE)) 
    res 
}

感謝您的提示來交叉驗證，它幫助了很多。

來源

2016-06-26 13:13:45 tester

您可以使用Fortran功能從

米勒，A. J.（1992年）減少計算時間頗多。算法AS 274：最小二乘例程爲補充了紳士。英國皇家統計學會社會。 C系列（應用統計），41（2），458-478。

您可以通過使用此代碼

# simulate data 
set.seed(101) 
n <- 1000 
X <- matrix(rnorm(10 * n), n, 10) 
y <- drop(10 + X %*% runif(10)) + rnorm(n) 
dat <- data.frame(y = y, X) 

# assign wrapper for biglm 
biglm_wrapper <- function(formula, data, min_window_size){ 
    mf <- model.frame(formula, data) 
    X <- model.matrix(terms(mf), mf) 
    y - model.response(mf) 

    n <- nrow(X) 
    p <- ncol(X) 
    storage.mode(X) <- "double" 
    storage.mode(y) <- "double" 
    w <- 1 
    qr <- list(
    d = numeric(p), rbar = numeric(choose(p, 2)), 
    thetab = numeric(p), sserr = 0, checked = FALSE, tol = numeric(p)) 
    nrbar = length(qr$rbar) 
    beta. <- numeric(p) 

    out <- numeric(n - min_window_size - 2) 
    for(i in 1:(n - 1)){ 
    row <- X[i, ] # will be over written 
    qr[c("d", "rbar", "thetab", "sserr")] <- .Fortran(
     "INCLUD", np = p, nrbar = nrbar, weight = w, xrow = row, yelem = y[i], 
     d = qr$d, rbar = qr$rbar, thetab = qr$thetab, sserr = qr$sserr, ier = i - 0L, 
     PACKAGE = "biglm")[ 
     c("d", "rbar", "thetab", "sserr")] 

    if(i >= min_window_size){ 
     coef. <- .Fortran(
     "REGCF", np = p, nrbar = nrbar, d = qr$d, rbar = qr$rbar, 
     thetab = qr$thetab, tol = qr$tol, beta = beta., nreq = p, ier = i, 
     PACKAGE = "biglm")[["beta"]] 
     out[i - min_window_size + 1] <- coef. %*% X[i + 1, ] 
    } 
    } 

    out 
} 

# assign function to compare with 
func <- function(formula, data, min_window_size){ 
    sapply(seq(min_window_size, nrow(data)-1), function(x) { 
    mod1 <- lm(formula, data = data[1:x, ]) 
    pred1 <- predict(mod1, newdata = data[x+1, ]) 
    pred1 
    }) 
} 

# show that the two gives the same 
frm <- y ~ X1 + X2 + X3 + X4 + X5 + X6 + X7 + X8 + X9 + X10 
r1 <- biglm_wrapper(frm, dat, 25) 
r2 <- func(frm, dat, 25) 
all.equal(r1, r2, check.attributes = FALSE) 
#R> [1] TRUE 

# show run time 
microbenchmark::microbenchmark(
    r1 = biglm_wrapper(frm, dat, 25), r2 = f2(frm, dat, 25), 
    times = 5) 
#R> Unit: milliseconds 
#R> expr   min   lq  mean  median   uq  max neval cld 
#R> r1 9.976505 10.00653 11.85052 10.53157 13.36377 15.37424  5 a 
#R> r2 1095.944410 1098.29661 1122.17101 1098.58264 1113.48833 1204.54306  5 b

來源

2018-02-18 14:14:25

創建遞歸式樣本外估計以有效計算R中的RMSE的方法

回答

相關問題