효율적으로 RMSE를 계산하기 위해 반복적 인 샘플 밖 추정을 만드는 방법

Question

다른 OLS 모델의 샘플 밖 예측 성능을 평가하려고합니다. 가장 쉬운 시계열 회귀는 다음과 같습니다 Y_t = B0 + B1 * Y_t-30 + e_t

이다가,의 50을 가정 해 봅시다 모델의 피팅 기간, 그때는 dynlm 패키지

를 사용하여 모델을 실행하자

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

현재 코드에서는 색인이 끝날 때까지 실행 한 다음 해당 모델의 RMSE를 저장합니다. 그러나이 RMSE는 샘플 밖의 한 단계 앞선 예측이 아니며 현재 코드가 이미 매우 느리고 원하는 작업을 정확하게 수행하지 못하기 때문에 사용자가 제안한 것이 있는지 물어보고 싶습니다. 패키지 내 목표를 달성하는 데 사용해야합니다. 앞서 한 단계를 만들 수)

1) 특정 피팅 기간이 지나면 재귀 회귀를 실행 (확대 창, NOT 롤 창)

2 :

, 나는 다음을 수행 할을 요약하면 아웃 오브 샘플은 루트 모델의 성능

지금까지와 함께이 일을 노력을 평가하기 위해 이러한 예측 대비 실제 관찰의 제곱 오차를 평균 계산)

3

예측 거대한에 대한 루프와 dynlm 패키지 ,하지만 res 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) 
} 

Answer 1

글쎄, 나는 내 문제 해결 방법에 기여하고 싶은 질문을 질문 한 것과 : 난 단지 내가 다른 모든 버릴 수있는 한 발 앞서 예측을 필요로

을이가 만든 코드를 빨리 실행하십시오. (모델 당 12 분 ~ 10 초).

저는 전체 데이터 프레임 (지연 포함)을 직접 생성하고 dynlm 대신 lm을 사용했습니다. 다음 코드는 내가 원하는 결과를 얻었습니다 (수동으로 처음 몇 관측을 확인했는데 작동하는 것처럼 보입니다). 코드는 여기에서 적응 : 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를 계산하는 내가 CrossValidated에 대한 힌트를

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

감사합니다, 그것은 많은 도움이 기능을 사용했다.

Answer 2

당신은

밀러, 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