搁肠辫辫のすすめ
?
29 likes?14,165 views
Rの関数をC++で記述することを可能にするパッケージRcppの使い方の紹介です。 第10回 Kashiwa.R で発表させて頂いた時の資料
![gibbsR <- function(N,thin){
mat<-matrix(0,nrow=N,ncol=2)
x <- 0
y <- 0
for(i in 1:N){
for(j in 1:thin){
x <- rgamma(1, 3, 1/(y*y+4))
y <- rnorm(1, 1/(x+1), 1/sqrt(2*x+2))
}
mat[i,] <- c(x,y)
}
return(mat)
}
Rのコード例
ギブスサンプラー
ついでに cmpfun でコンパイルしてみる
gibbsC <- compiler::cmpfun(gibbsR)](https://image.slidesharecdn.com/kashiwar10tsuda-140705090257-phpapp02/85/Rcpp-6-320.jpg)
![Rcppのコード例
#include <Rcpp.h>
using namespace Rcpp; ←これを書くとRcpp::は書かなくて良い
// [[Rcpp::export]] ←この直下の関数がRに読み込まれる
Rcpp::NumericMatrix gibbsCpp(int N, int thin) {
Rcpp::NumericMatrix mat(N, 2);//N行２列の行列を作成
double x = 0, y = 0;
for(int i = 0; i < N; ++i) {
for(int j = 0; j < thin; ++j) {
x = R::rgamma(3.0, 1.0 / (y * y + 4));
y = R::rnorm(1.0 / (x + 1), 1.0 / sqrt(2 * x + 2));
}
mat(i, 0) = x; // i 行 0 列に値を代入
mat(i, 1) = y; // i 行 1 列に値を代入
}
return(mat);//結果を返す
}
ギブスサンプラー](https://image.slidesharecdn.com/kashiwar10tsuda-140705090257-phpapp02/85/Rcpp-7-320.jpg)
![Vectorの作成と要素へのアクセス
要素番号でアクセス：v[0]?v(0)
最初の要素は0番、範囲外アクセスは落ちるので注意
要素の名前でアクセス：v[ 名前 ]
名前 が見つからない時は末尾に要素が追加される
NumericVector v(10); //0が10個
NumericVector::create(1,1,1); //1が３個
論理値ベクトルLでアクセス：v[L]
後で説明します
ベクターの作成](https://image.slidesharecdn.com/kashiwar10tsuda-140705090257-phpapp02/85/Rcpp-12-320.jpg)
![DataFrameの作成と要素へのアクセス
DataFrame::create( Named("列名1")= v1, _[ 列名2 ] = v2)
ベクトル v1, v2 からデータフレームを作成
v = df[0];
v = df[ 列名 ];
データフレーム df の列をベクター v に代入
ベクターに代入しないと
リストやデータフレームの値にはアクセス出来ない](https://image.slidesharecdn.com/kashiwar10tsuda-140705090257-phpapp02/85/Rcpp-14-320.jpg)
![Listの作成と要素へのアクセス
ベクターに代入しないと
リストやデータフレームの値にはアクセス出来ない
v = L[0];
v = L[ 要素名 ];
リストLの要素をベクターvに代入
ベクトル v1, v2 からリストを作成
List::create( Named("列名1")= v1, _[ 列名2 ] = v2)](https://image.slidesharecdn.com/kashiwar10tsuda-140705090257-phpapp02/85/Rcpp-15-320.jpg)
![Rcppオブジェクトの注意点
Rcpp
x<-1:5
myFunc(x)
x B A
A の変更は x に影響しない
NumericVector
myFunc(NumericVector A){
A[0]=100;//①
NumericVector B=A;//②
B[1]=200;
return(A);
}//A==c(100,200,3,4,5)B の変更が A に影響する?
変数→
①代入時に値がコピーされる
②
メモリー](https://image.slidesharecdn.com/kashiwar10tsuda-140705090257-phpapp02/85/Rcpp-16-320.jpg)
![Rcppオブジェクトの注意点
メモリー
NumericVector
myFunc(NumericVector A){
A[0]=100;//①
NumericVector B=clone(A);
B[1]=200;
return(A);
}//A==c(100,2,3,4,5)
x<-1:5
myFunc(x)
x
B の変更が A に影響しない
B A
clone()で値をコピーする
Rcpp
A の変更は x に影響しない
コピー
変数→
①代入時に値がコピーされる
②](https://image.slidesharecdn.com/kashiwar10tsuda-140705090257-phpapp02/85/Rcpp-17-320.jpg)
![?要素ごとにベクトル x と y の値を比較
?スカラーとベクトル要素の値を比較
LogicalVector res = (x==y);
x>y
x>=y
x!=y
論理値ベクトル
論理ベクトルでベクトルの要素にアクセスできる
??Vector z = x[res];
x[res] = val;
x==2
x!=2
x<2
x<=2](https://image.slidesharecdn.com/kashiwar10tsuda-140705090257-phpapp02/85/Rcpp-21-320.jpg)
![Rの関数を呼び出す
// [[Rcpp::export]]
NumericVector RcppFunc(Function f){
f(5, 1, 2);
//引数名を指定する場合
NumericVector res=f(5, Named("sd", 2),Named("mean", 1));
return(res);
}
１．Rcpp関数にRの関数を引数として渡す
２．Environmentクラス、Languageクラスを利用
//パッケージの環境を読み込む
Environment stats(“package:stats");
//関数を読み込む
Function rnorm = stats["rnorm"];
return rnorm(5, Named("sd", 2), Named("mean", 1));
Language call("rnorm", 5, Named("sd",2), Named("mean", 1));
return call.eval();](https://image.slidesharecdn.com/kashiwar10tsuda-140705090257-phpapp02/85/Rcpp-22-320.jpg)
![イテレータによる要素へのアクセス
NumericVector v;//合計を求める
double total = 0;
for(NumericVector::iterator i=v.begin(); i!=v.end(); ++i){
total += *i;
}
return(total);
メモリー
v[0] v[1] v[2] v[3]
i = v.begin() v.end()
++i
v[k]== *(v.begin()+k)
[ ]によるアクセスよりも一般に高速
++i ++i ++i](https://image.slidesharecdn.com/kashiwar10tsuda-140705090257-phpapp02/85/Rcpp-27-320.jpg)
![C++の標準ライブラリの利用
#include <numeric> //C++の標準ライブラリ
#include <Rcpp.h>
using namespace Rcpp;
!
// [[Rcpp::export]]
double Sum(NumericVector x) {//合計を求める
return std::accumulate(x.begin(), x.end(), 0.0);
}
C++には様々なアルゴリズムが用意されており
多くはイテレーターで位置や範囲を指定する
詳しくは cpprefjp などを参照
https://sites.google.com/site/cpprefjp/reference/algorithm](https://image.slidesharecdn.com/kashiwar10tsuda-140705090257-phpapp02/85/Rcpp-28-320.jpg)
![C++11を使う
Rcppコードの中で
// [[Rcpp::plugins("cpp11")]]
Sys.setenv( PKG_CXXFLAGS"="-std=c++11")
WindowsではRtoolsのgccバージョンがやや古いので
-std=c++0x と指定する
あるいはRで
auto、範囲for、イニシャライザリストなどを使うと
コーディングが楽になります](https://image.slidesharecdn.com/kashiwar10tsuda-140705090257-phpapp02/85/Rcpp-29-320.jpg)
![Boostを使う
Rで ?install.packages("BH")
Rcppで?// [[Rcpp::depends(BH)]]
その他の方法で Boost をインストールした場合
インストールしたパスを指定すればOK
BHパッケージにはBoostのうち
ヘッダーオンリーで使えるものが含まれている
Sys.setenv( PKG_CXXFLAGS =
"-I/opt/local/include -L/opt/local/lib/")](https://image.slidesharecdn.com/kashiwar10tsuda-140705090257-phpapp02/85/Rcpp-30-320.jpg)
![C++で演算子を実装する
NumericMatrix
operator+(NumericMatrix x, NumericMatrix y){
int nrow=x.nrow();
int ncol=x.ncol();
NumericMatrix res(nrow,ncol);
for(int i=0; i<nrow;++i){
for(int j=0; j<ncol;++j){
res(i,j)=x(i,j)+y(i,j);
}
}
return(res);
}
Rcppの記述の簡単さの秘密
=, [], (), >, == などあらゆる演算子を定義できる](https://image.slidesharecdn.com/kashiwar10tsuda-140705090257-phpapp02/85/Rcpp-31-320.jpg)
搁肠辫辫のすすめ
- 5. Rcpp使う C++コンパイラをインストール Windows：Rtools?Mac：Xcode?Linux：gcc RStudioで New File > C++ File Rで library(Rcpp) sourceCpp('ファイル名.cpp')
- 6. gibbsR <- function(N,thin){ mat<-matrix(0,nrow=N,ncol=2) x <- 0 y <- 0 for(i in 1:N){ for(j in 1:thin){ x <- rgamma(1, 3, 1/(y*y+4)) y <- rnorm(1, 1/(x+1), 1/sqrt(2*x+2)) } mat[i,] <- c(x,y) } return(mat) } Rのコード例 ギブスサンプラー ついでに cmpfun でコンパイルしてみる gibbsC <- compiler::cmpfun(gibbsR)
- 7. Rcppのコード例 #include <Rcpp.h> using namespace Rcpp; ←これを書くとRcpp::は書かなくて良い // [[Rcpp::export]] ←この直下の関数がRに読み込まれる Rcpp::NumericMatrix gibbsCpp(int N, int thin) { Rcpp::NumericMatrix mat(N, 2);//N行２列の行列を作成 double x = 0, y = 0; for(int i = 0; i < N; ++i) { for(int j = 0; j < thin; ++j) { x = R::rgamma(3.0, 1.0 / (y * y + 4)); y = R::rnorm(1.0 / (x + 1), 1.0 / sqrt(2 * x + 2)); } mat(i, 0) = x; // i 行 0 列に値を代入 mat(i, 1) = y; // i 行 1 列に値を代入 } return(mat);//結果を返す } ギブスサンプラー
- 8. library(rbenchmark) n <- 1000 thn <- 10 benchmark(gibbsR(n, thn), gibbsC(n, thn), gibbsCpp(n, thn), columns=c("test", "replications", "elapsed", "relative"), order="relative", replications=100) test replications elapsed relative 3 gibbsCpp(n, thn) 100 0.375 1.000 2 gibbsC(n, thn) 100 19.881 53.016 1 gibbsR(n, thn) 100 21.378 57.008 Rcpp版は50倍以上高速！！ 速度を比較
- 9. RとC++の違い ベクターなどの要素番号は 0 から 変数や関数（の返値）のデータ型を指定する必要がある 文の最後にセミコロン必須 変数名に . ドットは使えない、特別な意味がある 値の代入は =（-> は特別な意味がある） 関数の返値に return() 必須
- 10. Rcppの基本データ型 ??Matrix R Rcpp C++ logial Logical bool integer Integer int numeric Numeric double complex Complex complex character String string Date Date - POSIXct Datetime - 注：正確にはRcppに基本データ型として定義されているのは下の３つだけ 他は、ベクター、マトリックス型だけが定義されている
- 12. Vectorの作成と要素へのアクセス 要素番号でアクセス：v[0]?v(0) 最初の要素は0番、範囲外アクセスは落ちるので注意 要素の名前でアクセス：v[ 名前 ] 名前 が見つからない時は末尾に要素が追加される NumericVector v(10); //0が10個 NumericVector::create(1,1,1); //1が３個 論理値ベクトルLでアクセス：v[L] 後で説明します ベクターの作成
- 13. Matrixの作成と要素へのアクセス 行 ?列番号でアクセス：m( 1 , 2 ) ! NumericMatrix m(2, 2) 2行２列、値は0 i 行目をベクター v にコピー NumericVector v = m( i , _ ); 行?列の値を参照したいだけの場合 NumericMatrix::Row r=m( i , _ ); NumericMatrix::Column c=m( _ , j ); 作成 特定の行?列にアクセス：m( 1 , _ )
- 14. DataFrameの作成と要素へのアクセス DataFrame::create( Named("列名1")= v1, _[ 列名2 ] = v2) ベクトル v1, v2 からデータフレームを作成 v = df[0]; v = df[ 列名 ]; データフレーム df の列をベクター v に代入 ベクターに代入しないと リストやデータフレームの値にはアクセス出来ない
- 15. Listの作成と要素へのアクセス ベクターに代入しないと リストやデータフレームの値にはアクセス出来ない v = L[0]; v = L[ 要素名 ]; リストLの要素をベクターvに代入 ベクトル v1, v2 からリストを作成 List::create( Named("列名1")= v1, _[ 列名2 ] = v2)
- 16. Rcppオブジェクトの注意点 Rcpp x<-1:5 myFunc(x) x B A A の変更は x に影響しない NumericVector myFunc(NumericVector A){ A[0]=100;//① NumericVector B=A;//② B[1]=200; return(A); }//A==c(100,200,3,4,5)B の変更が A に影響する? 変数→ ①代入時に値がコピーされる ② メモリー
- 17. Rcppオブジェクトの注意点 メモリー NumericVector myFunc(NumericVector A){ A[0]=100;//① NumericVector B=clone(A); B[1]=200; return(A); }//A==c(100,2,3,4,5) x<-1:5 myFunc(x) x B の変更が A に影響しない B A clone()で値をコピーする Rcpp A の変更は x に影響しない コピー 変数→ ①代入時に値がコピーされる ②
- 18. ベクトルの長さなど属性値 ??Vector V V.names() V.length() ??Matrix M; M.ncol() M.nrow() M.attr( dimnames )= ??List::create(行名ベクタ, 列名ベクタ); DataFrame DF; DF.nrows() //行数 DF.size() //列数 DF.attr( names )//列名 DF.attr( row.names )//行名 List L; L.attr( names )//要素名
- 19. Rcppが提供するRライクな関数 ベクター関係 head, tail, rep_each, rep_len, rep, rev, seq_along, seq_len, cumsum, di?, pmin, pmax 確率分布 d/q/p/r 全てのRの提供する確率分布 数学関数 abs, acos, asin, atan, beta, ceil, ceiling, choose,cos, cosh, digamma, exp, expm1, factorial, ?oor, gamma, lbeta, lchoose, lfactorial, lgamma, log, log10, log1p, pentagamma, psigamma, round, signif, sin, sinh, sqrt, tan, tanh, tetragamma, trigamma, trunc 要約統計量 mean, min, max, sum, sd, var 値の検索 match, self_match, which_max, which_min 重複の値 duplicated, unique apply 関数 lapply, sapply, mapply 論理値 ifelse, all, any, is_true, is_false, is_na
- 20. Rcppのベクトル演算 NumericVector x, y; ベクトルの要素ごとの四則演算 x+y?x-y?x*y?x/y ベクトルとスカラーの乗除算 2*x?x/2 （v0.11.2現在）Matrix には対応してない ?????解決１：行列演算は RcppArmadillo ?????解決２：自分で演算子を実装する（後述）
- 21. ?要素ごとにベクトル x と y の値を比較 ?スカラーとベクトル要素の値を比較 LogicalVector res = (x==y); x>y x>=y x!=y 論理値ベクトル 論理ベクトルでベクトルの要素にアクセスできる ??Vector z = x[res]; x[res] = val; x==2 x!=2 x<2 x<=2
- 22. Rの関数を呼び出す // [[Rcpp::export]] NumericVector RcppFunc(Function f){ f(5, 1, 2); //引数名を指定する場合 NumericVector res=f(5, Named("sd", 2),Named("mean", 1)); return(res); } １．Rcpp関数にRの関数を引数として渡す ２．Environmentクラス、Languageクラスを利用 //パッケージの環境を読み込む Environment stats(“package:stats"); //関数を読み込む Function rnorm = stats["rnorm"]; return rnorm(5, Named("sd", 2), Named("mean", 1)); Language call("rnorm", 5, Named("sd",2), Named("mean", 1)); return call.eval();
- 23. Rcpp学習のためのサイト 本家のギャラリーには多くの例が掲載されている http://gallery.rcpp.org/ Advanced R の Rcpp のページはわかりやすいです http://adv-r.had.co.nz/Rcpp.html Rcpp Quick Reference Guide http://dirk.eddelbuettel.com/code/rcpp/Rcpp-quickref.pdf 手前味 ですが、Rcpp 入門（作成中） https://sites.google.com/site/rcppintroduction/ Rcppリファレンスはソースから自動生成されたもの http://dirk.eddelbuettel.com/code/rcpp/html/
- 24. C++学習のためのサイト http://ppp-lab.sakura.ne.jp/index_old.html Programming Place 沢山ありすぎるので 個人的にわかりやすかったものを少しだけ リファレンスは cpprefjp とか https://sites.google.com/site/cpprefjp/reference/algorithm
- 25. ?番外：颁++寄りな话题
- 26. RcppとC++データ構造の変換 vector<C++型> list<C++型> map<key, C++型> vector<vector<C++型>> List(Rcpp型Vector, …) Rcpp型Vector vector<double> Cpp = as(Rcpp); NumericVector Rcpp = wrap(Cpp); 変換可能なC++型はRcpp関数の引数や返値にもできる
- 27. イテレータによる要素へのアクセス NumericVector v;//合計を求める double total = 0; for(NumericVector::iterator i=v.begin(); i!=v.end(); ++i){ total += *i; } return(total); メモリー v[0] v[1] v[2] v[3] i = v.begin() v.end() ++i v[k]== *(v.begin()+k) [ ]によるアクセスよりも一般に高速 ++i ++i ++i
- 28. C++の標準ライブラリの利用 #include <numeric> //C++の標準ライブラリ #include <Rcpp.h> using namespace Rcpp; ! // [[Rcpp::export]] double Sum(NumericVector x) {//合計を求める return std::accumulate(x.begin(), x.end(), 0.0); } C++には様々なアルゴリズムが用意されており 多くはイテレーターで位置や範囲を指定する 詳しくは cpprefjp などを参照 https://sites.google.com/site/cpprefjp/reference/algorithm
- 29. C++11を使う Rcppコードの中で // [[Rcpp::plugins("cpp11")]] Sys.setenv( PKG_CXXFLAGS"="-std=c++11") WindowsではRtoolsのgccバージョンがやや古いので -std=c++0x と指定する あるいはRで auto、範囲for、イニシャライザリストなどを使うと コーディングが楽になります
- 30. Boostを使う Rで ?install.packages("BH") Rcppで?// [[Rcpp::depends(BH)]] その他の方法で Boost をインストールした場合 インストールしたパスを指定すればOK BHパッケージにはBoostのうち ヘッダーオンリーで使えるものが含まれている Sys.setenv( PKG_CXXFLAGS = "-I/opt/local/include -L/opt/local/lib/")
- 31. C++で演算子を実装する NumericMatrix operator+(NumericMatrix x, NumericMatrix y){ int nrow=x.nrow(); int ncol=x.ncol(); NumericMatrix res(nrow,ncol); for(int i=0; i<nrow;++i){ for(int j=0; j<ncol;++j){ res(i,j)=x(i,j)+y(i,j); } } return(res); } Rcppの記述の簡単さの秘密 =, [], (), >, == などあらゆる演算子を定義できる