आर

Question

मैं एक रोलिंग प्रतिगमन बहुत निम्न कोड के लिए इसी तरह चल रहा हूँ में एक रोलिंग खिड़की प्रतिगमन parallelize:

library(PerformanceAnalytics) 
library(quantmod) 
data(managers) 

FL <- as.formula(Next(HAM1)~HAM1+HAM2+HAM3+HAM4) 
MyRegression <- function(df,FL) { 
    df <- as.data.frame(df) 
    model <- lm(FL,data=df[1:30,]) 
    predict(model,newdata=df[31,]) 
} 

system.time(Result <- rollapply(managers, 31, FUN="MyRegression",FL, 
    by.column = FALSE, align = "right", na.pad = TRUE)) 

मैं कुछ अतिरिक्त प्रोसेसर मिल गया है, तो मैं parallelize के लिए एक रास्ता खोजने की कोशिश कर रहा हूँ रोलिंग विंडो। यह एक गैर रोलिंग प्रतिगमन था तो मैं आसानी से कार्यों का परिवार लागू का उपयोग कर इसे parallelize सकता है ...

Answer 1

स्पष्ट एक ताकि आप सूत्र आदि के प्रसंस्करण में सभी भूमि के ऊपर उठाना नहीं है lm() के बजाय lm.fit() उपयोग करने के लिए है ।

अद्यतन: तो जब मैंने कहा कि स्पष्ट मैं क्या कहना मतलब था blindingly स्पष्ट लेकिन भ्रामक लागू करना मुश्किल!

लगभग नगण्य का एक सा के बाद, मैं इस

library(PerformanceAnalytics) 
library(quantmod) 
data(managers) 

पहले चरण को एहसास है कि मॉडल मैट्रिक्स पहले से बनाए गए हो सकता है के साथ आया था, इसलिए हम ऐसा है और इसे वापस कन्वर्ट एक चिड़ियाघर वस्तु के लिए rollapply() साथ उपयोग करें:

+०१२३५१६४१०६:

mmat2 <- model.frame(Next(HAM1) ~ HAM1 + HAM2 + HAM3 + HAM4, data = managers, 
        na.action = na.pass) 
mmat2 <- cbind.data.frame(mmat2[,1], Intercept = 1, mmat2[,-1]) 
mmatZ <- as.zoo(mmat2) 

अब हम एक समारोह है कि प्रत्येक यात्रा पर डिजाइन मैट्रिक्स बनाए बिना lm.fit() रोजगार बड़े कार्य करने करने के लिए की आवश्यकता होगी

MyRegression2 <- function(Z) { 
    ## store value we want to predict for 
    pred <- Z[31, -1, drop = FALSE] 
    ## get rid of any rows with NA in training data 
    Z <- Z[1:30, ][!rowSums(is.na(Z[1:30,])) > 0, ] 
    ## Next() would lag and leave NA in row 30 for response 
    ## but we precomputed model matrix, so drop last row still in Z 
    Z <- Z[-nrow(Z),] 
    ## fit the model 
    fit <- lm.fit(Z[, -1, drop = FALSE], Z[,1]) 
    ## get things we need to predict, in case pivoting turned on in lm.fit 
    p <- fit$rank 
    p1 <- seq_len(p) 
    piv <- fit$qr$pivot[p1] 
    ## model coefficients 
    beta <- fit$coefficients 
    ## this gives the predicted value for row 31 of data passed in 
    drop(pred[, piv, drop = FALSE] %*% beta[piv]) 
} 

समय की तुलना:

> system.time(Result <- rollapply(managers, 31, FUN="MyRegression",FL, 
+         by.column = FALSE, align = "right", 
+         na.pad = TRUE)) 
    user system elapsed 
    0.925 0.002 1.020 
> 
> system.time(Result2 <- rollapply(mmatZ, 31, FUN = MyRegression2, 
+         by.column = FALSE, align = "right", 
+         na.pad = TRUE)) 
    user system elapsed 
    0.048 0.000 0.05 

कौन सा मूल के ऊपर एक बहुत ही उचित सुधार मिलता। और अब जांच कि जिसके परिणामस्वरूप वस्तुओं एक ही कर रहे हैं:

> all.equal(Result, Result2) 
[1] TRUE 

का आनंद लें!

Answer 2

आप एक अपघटन अपडेट करके रन के समय को कम कर सकते हैं। यह के बजाय प्रत्येक पुनरावृत्ति पर लागत उत्पन्न करता है जहां n आप विंडो चौड़ाई है। नीचे दो की तुलना करने के लिए एक कोड है। यह सी ++ में बहुत तेज़ी से कर रहा है लेकिन LINPACK dchud और dchdd आर के साथ शामिल नहीं हैं इसलिए आपको ऐसा करने के लिए एक पैकेज लिखना होगा। इसके अलावा, मैं पढ़ने याद है कि आप Linpack dchud और dchdd आर अद्यतन

library(SamplerCompare) # for LINPACK `chdd` and `chud` 
roll_forcast <- function(X, y, width){ 
    n <- nrow(X) 
    p <- ncol(X) 
    out <- rep(NA_real_, n) 

    is_first <- TRUE 
    i <- width 
    while(i < n){ 
    if(is_first){ 
     is_first <- FALSE 
     qr. <- qr(X[1:width, ]) 
     R <- qr.R(qr.) 

     # Use X^T for the rest 
     X <- t(X) 

     XtY <- drop(tcrossprod(y[1:width], X[, 1:width])) 
    } else { 
     x_new <- X[, i] 
     x_old <- X[, i - width] 

     # update R 
     R <- .Fortran(
     "dchud", R, p, p, x_new, 0., 0L, 0L, 
     0., 0., numeric(p), numeric(p), 
     PACKAGE = "SamplerCompare")[[1]] 

     # downdate R 
     R <- .Fortran(
     "dchdd", R, p, p, x_old, 0., 0L, 0L, 
     0., 0., numeric(p), numeric(p), integer(1), 
     PACKAGE = "SamplerCompare")[[1]] 

     # update XtY 
     XtY <- XtY + y[i] * x_new - y[i - width] * x_old 
    } 

    coef. <- .Internal(backsolve(R, XtY, p, TRUE, TRUE)) 
    coef. <- .Internal(backsolve(R, coef., p, TRUE, FALSE)) 

    i <- i + 1 
    out[i] <- X[, i] %*% coef. 
    } 

    out 
} 

# simulate data 
set.seed(101) 
n <- 10000 
wdth = 50 
X <- matrix(rnorm(10 * n), n, 10) 
y <- drop(X %*% runif(10)) + rnorm(n) 
Z <- cbind(y, X) 

# assign other function 
lm_version <- function(Z, width = wdth) { 
    pred <- Z[width + 1, -1, drop = FALSE] 
    ## fit the model 
    Z <- Z[-nrow(Z), ] 
    fit <- lm.fit(Z[, -1, drop = FALSE], Z[,1]) 
    ## get things we need to predict, in case pivoting turned on in lm.fit 
    p <- fit$rank 
    p1 <- seq_len(p) 
    piv <- fit$qr$pivot[p1] 
    ## model coefficients 
    beta <- fit$coefficients 
    ## this gives the predicted value for row 31 of data passed in 
    drop(pred[, piv, drop = FALSE] %*% beta[piv]) 
} 

# show that they yield the same 
library(zoo) 
all.equal(
    rollapply(Z, wdth + 1, FUN = lm_version, 
      by.column = FALSE, align = "right", fill = NA_real_), 
    roll_forcast(X, y, wdth)) 
#R> [1] TRUE 

# benchmark 
library(compiler) 
roll_forcast <- cmpfun(roll_forcast) 
lm_version <- cmpfun(lm_version) 
microbenchmark::microbenchmark(
    new = roll_forcast(X, y, wdth), 
    prev = rollapply(Z, wdth + 1, FUN = lm_version, 
        by.column = FALSE, align = "right", fill = NA_real_), 
    times = 10) 
#R> Unit: milliseconds 
#R> expr  min  lq  mean median  uq  max neval cld 
#R> new 113.7637 115.4498 129.6562 118.6540 122.4930 230.3414 10 a 
#R> prev 639.6499 674.1677 682.1996 678.6195 686.8816 763.8034 10 b