• こちらのマハラノビス距離の繰り返し計算でスピードアップの様子を見てみよう
• cxxfunction()関数を使っている例をsourceCpp()用に書き換えてみる
• データ型指定はRcpp namespace形式らしい
```## Using RcppArmadillo
code <- 'arma::mat x = Rcpp::as<arma::mat>(X);
arma::mat mu = Rcpp::as<arma::mat>(Mu);
arma::mat sigma = Rcpp::as<arma::mat>(Sigma);
int n = x.n_rows;
arma::vec md(n);
for (int i=0; i<n; i++){
arma::mat x_i = x.row(i) - mu;
arma::mat Y = arma::solve( sigma, arma::trans(x_i) );
md(i) = arma::as_scalar(x_i * Y);
}
return wrap(md);'
mahalanobis_RcppArma <- cxxfunction( signature(X="numeric",
```
```#include <RcppArmadilloExtensions/sample.h>

using namespace Rcpp ;

// [[Rcpp::export]]

NumericVector mahalanobis_RcppArma(NumericMatrix X,NumericMatrix Mu, NumericMatrix Sigma){
arma::mat x = Rcpp::as<arma::mat>(X);
arma::mat mu = Rcpp::as<arma::mat>(Mu);
arma::mat sigma = Rcpp::as<arma::mat>(Sigma);
int n = x.n_rows;
arma::vec md(n);
for (int i=0; i<n; i++){
arma::mat x_i = x.row(i) - mu;
arma::mat Y = arma::solve( sigma, arma::trans(x_i) );
md(i) = arma::as_scalar(x_i * Y);
}
return wrap(md);
}
```
```## Using R
mahalanobis_R <- function(X, Mu, Sigma){
md <- rep(NA, nrow(X) )
for(i in 1:nrow(X) )
md[i] <- mahalanobis(X[i,], center=Mu, cov=Sigma)
return(md)
}
```
```X <- matrix(runif(10000),100,100)
Mu <- matrix(runif(100),1,100)
Sigma <- matrix(runif(10000),100,100)

benchmark(
mahalanobis_R(X, c(Mu), Sigma),
mahalanobis_RcppArma(X, Mu, Sigma),
columns=c("test", "replications","elapsed", "relative"),
replications=100,
order="relative")
```
```                                test replications elapsed relative
2 mahalanobis_RcppArma(X, Mu, Sigma)          100   10.11    1.000
1     mahalanobis_R(X, c(Mu), Sigma)          100   49.02    4.849
```