Возможно, это поможет вам сократить время, значительно fastLm () намного медленнее, чем lm () ;Немного измените код fLmSEXP
, чтобы иметь возможность извлекать остатки.
library(Rcpp)
library(RcppArmadillo)
library(rbenchmark)
## start from SEXP, most conversions, longest code
src <- '
Rcpp::List fLmSEXP(SEXP Xs, SEXP ys) {
Rcpp::NumericMatrix Xr(Xs);
Rcpp::NumericVector yr(ys);
int n = Xr.nrow(), k = Xr.ncol();
arma::mat X(Xr.begin(), n, k, false);
arma::colvec y(yr.begin(), yr.size(), false);
// fit model y ~ X, extract residuals
arma::colvec coef = arma::solve(X, y);
arma::colvec res = y - X*coef;
// return the results
return Rcpp::List::create(Rcpp::Named("coefficients")=coef,Rcpp::Named("res")=res);
}
'
cppFunction(code=src, depends="RcppArmadillo")
Я создаю свой фрейм данных
df <- data.frame(replicate(3,sample(1:4,300000,rep=TRUE)))
df = cbind(X = rnorm(300000),df)
head(df)
X X1 X2 X3
1 0.6269854 1 4 3
2 0.4641201 1 1 4
3 -0.5625020 3 1 4
4 0.0452215 2 1 2
5 2.2453335 3 3 2
6 0.4045328 1 3 3
m <- as.matrix(cbind(X = df[,1],cbind(I = 1,df[,2:4])))
Я сравниваю результаты обеих функций
benchmark(
lm_res = lm(X ~ X1 + X2 + X3, data = df)$residuals,
flm_res = fLmSEXP(m[,2:5],m[,1])$res, replications = 100)[,1:4]
test replications elapsed relative
2 flm_res 100 4.14 1.00
1 lm_res 100 12.46 3.01
Надеюсь, это будет полезно или, по крайней мере, поможет вам.