Armadilloを使う〜RstudioでC++を覚える - ryamadaの遺伝学・遺伝統計学メモ

Rcppパッケージの基本、それをRstudioで扱うことの基本はわかった
けれど、「素のC++」だけでは、ごつすぎて使いにくい。特に線形代数とか、いわゆる数値計算のためのベクトル・行列・アレイなどが使いやすくて、それの演算ライブラリがないとだめなので、Armadilloを入れて、それを使えるようになることにする
そもそもArmadilloはC++用の線形代数ライブラリ(こちら)
それをRcppで使うようにしたのがRcppArmadilloパッケージなので、以下をして使えるようにする

install.packages("RcppArmadillo")
library(RcppArmadillo)

RcppArmadilloのチュートリアルとかもあるけれど、手っ取り早く、Rstudioで動かしているサンプルを見つけよう→こちら
- 結局、以下を"LinearRegression.cpp"という名前で保存して(これより上にヘッダーは不要(#include
  も不要)
  - 逆行列計算ができていたり(solve)、行列の転置ができていたり(mX.t())、行列の積ができていたり(mX.t() * mX)するところがArmadilloライブラリの恩恵。List型でRに返したりしているところはRcppArmadilloパッケージの実装の恩恵

// LinearRegression.cpp
// [[Rcpp::depends(RcppArmadillo)]]
#include <RcppArmadillo.h>
using namespace arma;  // use the Armadillo library for matrix computations
using namespace Rcpp;

// [[Rcpp::export]]
List fnLinRegRcpp(vec vY, mat mX) {

  // compute the OLS estimator & model residuals
  vec vBeta =  solve(mX.t()*mX, mX.t()*vY);  
  vec vResid = vY - mX * vBeta;

  // construct the return object
  List ret;
  ret["beta"] = vBeta;
  ret["resid"] = vResid;

  return ret;
}

// END

- R側でこのファイルを読み込んでコンパイルして("sourceCpp("hoge.cpp")"の代りに"Rcpp::sourceCpp("hoge.cpp")"

library(devtools)
library(Rcpp)
library(RcppArmadillo)
Rcpp::sourceCpp("LinearRegression.cpp", showOutput = TRUE, rebuild = FALSE) # もちろんファイルの置き場所はこれでアクセスできるように…

- 関数を使う

# generate some sample data
iK = 4
iN = 100
mX = cbind(1, matrix(rnorm(iK*iN), iN, iK))
vBeta0 = c(2, 3.5, 0.11, 6.33, 23)
vY = rnorm(iN, mean = mX %*% vBeta0)

# test the function
linReg1 = fnLinRegRcpp(vY, mX)
linReg1$beta  # coefficient estimates

# compare the results to the built-in lm.fit function
lm.fit(y = vY, x = mX)$coef  # coefficient estimates

# END