c++ - 为什么RcppArmadillo的fastLmPure在输出中产生NA，而fastLm没有？

我在R中大量使用了滚动回归，而我的初始设置是这样的:

dolm <- function(x) coef(lm(x[,1] ~ x[,2] + 0, data = as.data.frame(x)))
rollingCoef = rollapply(someData, 100, dolm)

上面的示例工作得很好，但是如果您有很多迭代的话它会很慢。

为了加快速度，我决定尝试使用Rcpp包。

首先，我用lm替换了fastLm，结果有点快，但仍然很慢。因此，这促使我尝试在c++中为循环编写整个滚动回归的系数函数，而不是在Rcpp帮助下将其集成到R中。

因此，我已将原始RcppArmadillo的函数fastLm更改为此:

// [[Rcpp::depends(RcppArmadillo)]]

#include <RcppArmadillo.h>

using namespace Rcpp;

// [[Rcpp::export]]

List rollCoef(const arma::mat& X, const arma::colvec& y, double window ) {

    double cppWindow = window - 1;
    double matRows = X.n_rows;
    double matCols = X.n_cols - 1;

    arma::mat coef( matRows - cppWindow, X.n_cols);   // matrix for estimated coefficients

    //for loop for rolling regression.
    for( double i = 0 ; i < matRows - cppWindow ; i++  )
    {
        coef.row(i) = arma::trans(arma::solve(X( arma::span(i,i + cppWindow), arma::span(0,matCols)) , y.rows(i,i + cppWindow)));
    }

  return List::create(_["coefficients"] = coef);
}

然后用sourceCpp(file=".../rollCoef.cpp")将其下载到R

因此，它比rollapply快得多，并且在较小的示例上也能很好地工作，但是比我将它应用于〜200000观察到的数据时，它在输出中产生了大约一半的NA，而rollapply / fastLm组合却什么也没有。

所以在这里我需要一些帮助。我的功能出了什么问题？为什么我的函数输出中没有NA，而rollapply / fastLm中没有NA，但是，如果我理解正确，它们都基于arma::solve？非常感谢您的帮助。

更新
这是可复制的代码:

require(Rcpp)
require(RcppArmadillo)
require(zoo)
require(repmis)
myData <- source_DropboxData(file = "example.csv",
                              key = "cbrmkkbssu5bn96", sep = ",", header = TRUE)

## in order to use my custom function "rollCoef" you should download it to R.
## The c++ code is presented above in the main question.
## Download it where you want as "rollCoef.cpp" and then download it to R with:

sourceCpp(file=".../rollCoeff.cpp"). # there should be your actual path.

myCoef = rollCoef(as.matrix(myData[,2]),myData[,1],260)

summary(unlist(myCoef)) # 80923 NA's

dolm = function(x) coef(fastLmPure(as.matrix(x[,2]), x[,1]))

myCoef2 = rollapply(myData, 260, dolm, by.column = FALSE)

summary(myCoef2) # 80923 NA's

dolm2 = function(x) coef(fastLm(x[,1] ~ x[,2] + 0, data = as.data.frame(x)))

myCoef3 = rollapply(myData, 260, dolm2, by.column = FALSE)

summary(myCoef3) # !!! No NA's !!!

head(unlist(myCoef)) ; head(unlist(myCoef2)) ; head(myCoef3)

因此，我的函数的输出与RcppArmadillo的fastLmPure和rollapply组合的输出相同，并且它们均产生NA，但带有rollapply的fastLm却不。据我了解，例如从HERE和HERE，fastLm基本上是在调用fastLmPure，但是为什么在第三种方法中没有NA呢？ fastLm中是否有一些其他功能可以阻止我没有发现的NA？

最佳答案

有一个完整的RcppRoll包可以进行自定义滚动-您应该能够对其进行扩展以及它的rollit()函数也可以滚动lm()。