1
2library(Rcpp)
3
4cppFunction('
5    NumericVector convolveCpp(NumericVector a, NumericVector b) {
6
7        int na = a.size(), nb = b.size();
8        int nab = na + nb - 1;
9        NumericVector xab(nab);
10
11        for (int i = 0; i < na; i++)
12            for (int j = 0; j < nb; j++)
13                xab[i + j] += a[i] * b[j];
14
15        return xab;
16    }
17')
18
19convolveCpp(c(1,2,3), matrix(3,3))
20
21
22cppFunction(depends='RcppArmadillo', code='
23    List fastLm(NumericVector yr, NumericMatrix Xr) {
24
25        int n = Xr.nrow(), k = Xr.ncol();
26
27        arma::mat X(Xr.begin(), n, k, false); // reuses memory and avoids copy
28        arma::colvec y(yr.begin(), yr.size(), false);
29
30        arma::colvec coef = arma::solve(X, y);      // fit model y ~ X
31        arma::colvec resid = y - X*coef;            // residuals
32
33        double sig2 = arma::as_scalar( arma::trans(resid)*resid/(n-k) );
34                                                    // std.error of estimate
35        arma::colvec stderrest = arma::sqrt(
36                        sig2 * arma::diagvec( arma::inv(arma::trans(X)*X)) );
37
38        return List::create(Named("coefficients") = coef,
39                            Named("stderr")       = stderrest
40    );
41}
42')
43
44