Почему тест Портманто медленнее в Rcpp, чем в R?

Question

Мне нужно провести тщательное изучение основных тестов Portmanteau Статья , для этого я должен оценить их в различных сценариях, размерах выборки и в разных моделях ARMA (p, q), генерируя 180 сценариев, что приводит меня кзакрыть 6 часов.Запрограммируйте мою функцию в R и Rcpp, однако я нахожу удивительным, что в C ++ она медленнее, мой вопрос почему?

My R Code:

Portmanteau <- function(x,h=1,type = c("Box-Pierce","Ljun-Box","Monti"),fitdf = 0){
  Ti <- length(x)
  df <- h-fitdf
  ri <- acf(x, lag.max = h, plot = FALSE, na.action = na.pass)
  pi <- pacf(x, lag.max = h, plot = FALSE, na.action = na.pass)
  if(type == "Monti"){d<-0} else{d<-1}
  if(type == "Box-Pierce"){wi <- 1} else{wi <- (Ti+2)/seq(Ti-1,Ti-h)}
  Q <- Ti*(d*sum(wi*identity(ri$acf[-1]^2))+(1-d)*sum(wi*identity(pi$acf^2)))
  pv <- pchisq(Q,df,lower.tail = F)
  result <- cbind(Statistic = Q, df,p.value = pv)
  rownames(result) <- paste(type,"test")
  return(result)
  }

My Rcpp code

#include <Rcpp.h>
using namespace Rcpp;

// [[Rcpp::export]]
NumericVector PortmanteauC(NumericVector x, int h = 1,const char* type = "Box-Pierce" ,int fitdf = 0) {
  Environment stats("package:stats");
  Function acf = stats["acf"];
  Function pacf = stats["pacf"];
  Function na_pass = stats["na.pass"];
  List ri =  acf(x, h, "correlation", false, na_pass);
  List pi =  pacf(x, h, false, na_pass);
  int Ti = x.size();
  int df = h - fitdf;
  double d; 
  NumericVector wi;
  NumericVector rk = ri["acf"];
  NumericVector pk = pi["acf"];
  NumericVector S(h);
  for(int i = 0; i < h; ++i){S[i] = Ti-i-1;}
  rk.erase(0);
  if(strcmp(type,"Monti") == 0){d=0;} else{d=1;}
  if(strcmp(type,"Box-Pierce") == 0){wi = rep(1,h);} else{wi = (Ti+2)/S;}
  double Q = Ti*(d*sum(wi*pow(rk,2)) + (1-d)*sum(wi*pow(pk,2)));
  double pv = R::pchisq(Q,df,0,false);
  NumericVector result(3);
  result[0] = Q;
  result[1] = df;
  result[2] = pv;
  return(result);
}

Пример

set.seed(1)
y = arima.sim(model = list(ar = 0.5), n = 250)
mod = arima(y, order = c(1,0,0))
res = mod$residuals

---

Box-Pierce

library(rbenchmark)
benchmark(PortmanteauC(res, h=10, type = "Box-Pierce",fitdf = 1),replications = 500,Portmanteau(res,h = 10, type = "Box-Pierce", fitdf= 1),
    Box.test(res, lag = 10, type = "Box-Pierce", fitdf= 1))[,1:4]

                                                       test replications elapsed relative
3   Box.test(res, lag = 10, type = "Box-Pierce", fitdf = 1)          500    0.17    1.000
2  Portmanteau(res, h = 10, type = "Box-Pierce", fitdf = 1)          500    0.44    2.588
1 PortmanteauC(res, h = 10, type = "Box-Pierce", fitdf = 1)          500    1.82   10.706

---

Ljun-Box

benchmark(Box.test(res, lag = 5, type = "Ljung-Box", fitdf= 1),replications = 500,
Portmanteau(res,h = 5, type = "Ljung-Box", fitdf= 1),
PortmanteauC(res,h = 5, type = "Ljung-Box", fitdf= 1))[,1:4]
                                                     test replications elapsed relative
1   Box.test(res, lag = 5, type = "Ljung-Box", fitdf = 1)          500    0.17    1.000
2  Portmanteau(res, h = 5, type = "Ljung-Box", fitdf = 1)          500    0.45    2.647
3 PortmanteauC(res, h = 5, type = "Ljung-Box", fitdf = 1)          500    1.84   10.824

---

Я бы ожидал, что Rcpp будет намного быстрее, чем скомпилированный R.

Answer 1

Давайте проанализируем свойства производительности вашего R-кода.Поскольку индивидуальный вызов настолько быстр, что профилировщик выборки, предоставленный R, не может быть легко использован, я просто использую repeat(), чтобы повторять код до тех пор, пока он не прервется:

Portmanteau <- function(x,h=1,type = c("Box-Pierce","Ljun-Box","Monti"),fitdf = 0){
  Ti <- length(x)
  df <- h-fitdf
  ri <- acf(x, lag.max = h, plot = FALSE, na.action = na.pass)
  pi <- pacf(x, lag.max = h, plot = FALSE, na.action = na.pass)
  if(type == "Monti"){d<-0} else{d<-1}
  if(type == "Box-Pierce"){wi <- 1} else{wi <- (Ti+2)/seq(Ti-1,Ti-h)}
  Q <- Ti*(d*sum(wi*identity(ri$acf[-1]^2))+(1-d)*sum(wi*identity(pi$acf^2)))
  pv <- pchisq(Q,df,lower.tail = F)
  result <- cbind(Statistic = Q, df,p.value = pv)
  rownames(result) <- paste(type,"test")
  return(result)
}

set.seed(1)
profvis::profvis({
  repeat({
    y = arima.sim(model = list(ar = 0.5), n = 250)
    mod = arima(y, order = c(1,0,0))
    res = mod$residuals
    Portmanteau(res, h = 10, type = "Box-Pierce", fitdf = 1)
  })
})

Я позволил ему работать в течение примерно 49 с.Часть графического вывода, предоставленного в RStudio, можно увидеть здесь:

Мы учимся на этом:

• arima() занимает в семь раз больше времени, чем Portmenteau().В зависимости от соотношения вызовов между этими двумя функциями вы можете оптимизировать неправильную функцию.
• Для вызова Portmenteau() почти полное время затрачивается на pacf() и acf().Эти функции R также используются в вашем коде Rcpp, но с дополнительными сложностями возврата к R из C ++.Это объясняет, почему ваш C ++ медленнее, чем ваш код R.