Rcpp, RcppArmadilloの練習(1)

## 0. Rcpp,RcppArmadilloの環境構築

実行環境はMac．MacでRcpp,RcppArmadilloをインストールするのはだいぶめんどくさい．

1. Rでinstall.packages("Rcpp", dependencies=TRUE)
2. XcodeとCommand Line Toolをインストールしてから，ターミナルで「sudo xcodebuild -license」と打ち込んでライセンス条項にagreeする．XcodeのインストールはApp Storeから．Command Line Toolは[こちら](http://hiroki.jp/mavericks-command-line-developer-tools)を参照．
3. Rでinstall.packages("RcppArmadillo", type = "source")
4. 3.で「library not found for -lgfortran
clang: error: linker command failed with exit code 1 (use -v to see invocation)」みたいなエラーが出る．gfortanが入っていない場合は[こちら](https://gcc.gnu.org/wiki/GFortranBinariesMacOS)からインストールする
5. 4.のエラーが出るのにgfortanが入っているはずの場合は，Finderからメニューの「移動」→「フォルダへ移動」と選んで，「~/.R」と検索する．ここに「Makevars」というファイルがあればそれを開き，無ければテキストエディタで新規作成する．Makevarsファイルに入れたら，「FLIBS=-L/usr/local/gfortran」という行を加える．
6. もう一度install.packages("RcppArmadillo", type = "source")を試し，エラーが出なければ完了．


## 1.行列の基本演算 

以下，[こちら](https://sites.google.com/view/omori-stat/japanese/ソフトウェア/rcpp?authuser=0)を参考にまとめる．

```
// MySample1.cpp
#include <RcppArmadillo.h> // 呪文
// [[Rcpp::depends(RcppArmadillo)]] // 呪文
// [[Rcpp::plugins("cpp11")]] // 呪文
// [[Rcpp::export]] // 呪文
void fSample1() // 関数fSample1を定義。戻り値なし。
{
  
  double da, db, dc; // 実数 
  arma::mat mA, mB, mC, mD; // RcppArmadilloに依存する行列
  
  
  /*---実数の定義---*/
  da = 2.0;
  db = arma::datum::pi; // piとかeとか特別な数をdatumで定義
  dc = da + db;
  
  Rprintf("da=%.4f, db=%.4f\n",da,db); // 表示にはRprintfかRcoutを使う。有効数字4桁
  Rprintf("da+db=%.4f\n",da+db);
  Rcpp::Rcout << "dc=" << dc <<"\n"; // オブジェクトを<<でつなげる
  
  
  /*---行列の定義---*/
  mA = { {1, 2}, 
         {3, 4} };
  
  mA.print("mA="); // 行列の表示
  
  mB = arma::zeros(2,2); //ゼロ行列 
  
  mB(0,0) = 10;
  mB(0,1) = 20;
  mB(1,0) = 30;
  mB(1,1) = 40;
  mB.print("mB=");
  
  mC ={ {100,200}, 
        {300,400}};
  mC.print("mC=");
  
  
  /*---簡単な行列演算---*/
  (da*mA).print("a*mA=");
  (mA+mB).print("mA+mB=");
  (mA*mB).print("mAmB="); // 行列積

  mD = arma::inv(mA); // 逆行列
  mD.print("mA^{-1}=");
  
  
  /*---ベクトルの定義---*/
  arma::vec vx = {1, 4, 3, 2, 5};
  vx.print("vx=");
  
  
  /*---簡単なベクトル演算---*/
  Rprintf("Mean     = %.4f \n", arma::mean(vx));
  Rprintf("Variance = %.4f \n", arma::var(vx));
  Rprintf("Minimum  = %.4f \n", arma::min(vx));
  Rprintf("Maximum  = %.4f \n", arma::max(vx));
  
  sort(vx, "ascend").print("Sort in ascending order");
  sort(vx, "descend").print("Sort in descending order");

  
  /*---簡単な行列演算part2---*/
  arma::mat mX = { {1, 30}, {4, 20}, {3, 50}, {2, 40}, {5, 10} };
  mX.print("mX=");
  arma::mean(mX).print("Mean");
  arma::var(mX).print("Variance");
  arma::min(mX).print("Minimum");
  arma::max(mX).print("Maximum");
  
  sort(mX, "descend").print("Sort in descending order");
}

/*** R
library(Rcpp)
library(RcppArmadillo)
fSample1() 
*/
```

```
###出力###
da=2.0000, db=3.1416
da+db=5.1416
dc=5.14159
mA=
   1.0000   2.0000
   3.0000   4.0000
mB=
   10.0000   20.0000
   30.0000   40.0000
mC=
   1.0000e+02   2.0000e+02
   3.0000e+02   4.0000e+02
a*mA=
   2.0000   4.0000
   6.0000   8.0000
mA+mB=
   11.0000   22.0000
   33.0000   44.0000
mAmB=
   7.0000e+01   1.0000e+02
   1.5000e+02   2.2000e+02
mA^{-1}=
  -2.0000   1.0000
   1.5000  -0.5000
vx=
   1.0000
   4.0000
   3.0000
   2.0000
   5.0000
Mean     =   3.0000 
Variance =   2.5000 
Minimum  =   1.0000 
Maximum  =   5.0000 
Sort in ascending order
   1.0000
   2.0000
   3.0000
   4.0000
   5.0000
Sort in descending order
   5.0000
   4.0000
   3.0000
   2.0000
   1.0000
mX=
    1.0000   30.0000
    4.0000   20.0000
    3.0000   50.0000
    2.0000   40.0000
    5.0000   10.0000
Mean
    3.0000   30.0000
Variance
   2.5000e+00   2.5000e+02
Minimum
    1.0000   10.0000
Maximum
    5.0000   50.0000
Sort in descending order
    5.0000   50.0000
    4.0000   40.0000
    3.0000   30.0000
    2.0000   20.0000
    1.0000   10.0000
```

## 2. 簡単な行列の作り方

```
// MySample2.cpp
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// [[Rcpp::plugins("cpp11")]]
// [[Rcpp::export]]
void fSample2() // 関数fSample2を定義，戻り値なし
{
  arma::mat mX, mY;
  arma::vec va, vb;
  
  
  /*---ゼロ行列，単位行列の定義---*/
  mX = arma::zeros(2,2);
  mY = arma::eye(2,2);
  
  mX.print("mX=");
  mY.print("mY=");
  
  
  /*---特別なベクトルの定義---*/
  va = arma::ones(2); // 1のみを要素とする長さ2のベクトル
  vb = arma::regspace(1,0.5,5); // 1から5まで,0.5刻みの等差数列
  
  va.print("va=");
  vb.print("vb=");
  (mY*va).print("mY*va=");
  
  
  /*---行列同士の操作---*/
  arma::vec vp, vq, vr;

  vp = {1, 2};
  vq = arma::zeros(2);
  vq[0] = 3;
  vq[1] = 4;

  vp.print("vp=");
  vq.print("vq=");
  
  arma::mat mP, mQ;
  
  mP = join_horiz(vp, vq); // 行方向に連結(2行2列になる)
  mQ = join_vert(vp, vq);  // 列方向に連結(4行1列になる)
  
  mP.print("mP=(vp vq)=");
  mQ.print("mQ=(vp' vq')'=");
  
  (vp.t()*vq).print("vp'vq="); // ベクトル積を計算.vp.t()で転置．．
  (arma::as_scalar(vp.t()*vq)*mP).print("(vp'*vq)*mP="); // as.scalarで一次元ベクトルを数値化
  
  }

/*** R
 fSample2()
*/
  
```

```
mX=
        0        0
        0        0
mY=
   1.0000        0
        0   1.0000
va=
   1.0000
   1.0000
vb=
   1.0000
   1.5000
   2.0000
   2.5000
   3.0000
   3.5000
   4.0000
   4.5000
   5.0000
mY*va=
   1.0000
   1.0000
vp=
   1.0000
   2.0000
vq=
   3.0000
   4.0000
mP=(vp vq)=
   1.0000   3.0000
   2.0000   4.0000
mQ=(vp' vq')'=
   1.0000
   2.0000
   3.0000
   4.0000
vp'vq=
   11.0000
(vp'*vq)*mP=
   11.0000   33.0000
   22.0000   44.0000
```

## 3.行列の要素を取り出す

```
// MySample3.cpp
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// using namespace arma;
// [[Rcpp::export]]

void fSample3()
{
  arma::mat mA, mB;
  double da;

  
  /*---行列の要素を参照---*/
  mA = { {1,2,3}, 
         {4,5,6},
         {7,8,9}};
  mA.print("mA=");
  ( mA.t() ).print("mA'=");
  
  da = mA(0, 0); // 行列mAの[0,0]要素をdaとおく．
  Rprintf("mA[0][0]= %5.3f\n", da);
  Rprintf("mA[0][1]= %5.3f\n", mA(0,1));
  Rprintf("mA[1][1]= %5.3f\n", mA(1,1));

  
  /*---スライス---*/
  ( mA(arma::span(0,0), arma::span(0,2)) ).print("mA[0][0:2]");
  ( mA(arma::span(0,2), arma::span(0,0)) ).print("mA[0:2][0]");
  ( mA(arma::span(0,1), arma::span(1,2)) ).print("mA[0:1][1:2]");
  ( mA(arma::span(1,2), arma::span(1,2)) ).print("mA[1:2][1:2]");
  
}  


/*** R
fSample3()
*/
```

```
###出力###
mA=
   1.0000   2.0000   3.0000
   4.0000   5.0000   6.0000
   7.0000   8.0000   9.0000
mA'=
   1.0000   4.0000   7.0000
   2.0000   5.0000   8.0000
   3.0000   6.0000   9.0000
mA[0][0]= 1.000
mA[0][1]= 2.000
mA[1][1]= 5.000
mA[0][0:2]
   1.0000   2.0000   3.0000
mA[0:2][0]
   1.0000
   4.0000
   7.0000
mA[0:1][1:2]
   2.0000   3.0000
   5.0000   6.0000
mA[1:2][1:2]
   5.0000   6.0000
   8.0000   9.0000
```

## 行列のインクリーメント，クロネッカー積

```
// MySample4.cpp
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// using namespace Rcpp
// using namespace arma;
// [[Rcpp::plugins("cpp11")]]
// [[Rcpp::export]]
void fSample4a()
{
 arma:: mat mA, mB, vc;
  
  mA = {{1,2},
        {3,4}};
  mA.print("mA=");
  mB = mA; // 代入．mAをいじってもmBは変わらない．

  mA++;  // mA = mA + 1; or mA += 1;と同値．各要素に1を加える．
  mA.print("mA=");
  mB.print("mB=");
  
  vc = {10, 20}; // 1 × 2 vector (縦ベクトル)
  vc = vc.t();  // 2 x 1 vector (横ベクトル)
  vc.print("vc=");
  
  mA--;  // mA = mA - 1; or mA -= 1;と同値．各要素から1引く．
  mA.print("mA=");
  
  (mA*mA).print("mA^2="); // 行列積
  (mA%mA).print("mA.^2="); // 要素積
  
  mB++;
  mB.print("mB=");
  ( mA*arma::inv(mB) ).print("mA*mB^{-1}="); // 行列商
  ( mA/mB).print("mA./mB="); // 各要素の割り算
  
  (mA*mB).print("mA*mB="); // 行列積
  (mA%mB).print("mA.*mB="); // 要素積
  ( arma::kron(mA, vc.t()) ).print("mA**vc'"); // クロネッカー積
}

// [[Rcpp::export]]
void fSample4b()
{
  double dx1, dx2, dx3, dx4;
  dx1 = 0; dx2= 1; dx3 = 2; dx4= 3;
  Rprintf("x1=%5.3f, x2=%5.3f, x3=%5.3f, x4=%5.3f\n", dx1, dx2, dx3, dx4);
  
  dx1 +=1; dx2 -=1; dx3 *=2, dx4 /=2;
  Rprintf("x1=%5.3f, x2=%5.3f, x3=%5.3f, x4=%5.3f\n", dx1, dx2, dx3, dx4);
}

/*** R
fSample4a() 
fSample4b() 
*/

```

```
###出力###
mA=
   1.0000   2.0000
   3.0000   4.0000
mA=
   2.0000   3.0000
   4.0000   5.0000
mB=
   1.0000   2.0000
   3.0000   4.0000
vc=
   10.0000
   20.0000
mA=
   1.0000   2.0000
   3.0000   4.0000
mA^2=
    7.0000   10.0000
   15.0000   22.0000
mA.^2=
    1.0000    4.0000
    9.0000   16.0000
mB=
   2.0000   3.0000
   4.0000   5.0000
mA*mB^{-1}=
   1.5000  -0.5000
   0.5000   0.5000
mA./mB=
   0.5000   0.6667
   0.7500   0.8000
mA*mB=
   10.0000   13.0000
   22.0000   29.0000
mA.*mB=
    2.0000    6.0000
   12.0000   20.0000
mA**vc'
   10.0000   20.0000   20.0000   40.0000
   30.0000   60.0000   40.0000   80.0000

> fSample4b() 
x1=0.000, x2=1.000, x3=2.000, x4=3.000
x1=1.000, x2=0.000, x3=4.000, x4=1.500
```


## 5. Loop

```
// MySample5.cpp
#include <RcppArmadillo.h>
// [[Rcpp::depends(RcppArmadillo)]]
// using namespace arma;

// [[Rcpp::export]]
void fSample5()
{
  int i;
  
  
  /*---for文---*/
  Rprintf("loop: increasing\n");
  for( i=0; i < 3; ++i)
    {
      Rprintf("i=%i\n",i);
    }
  
  Rprintf("loop: decreasing\n");
  for( i=3; i > 0; --i)
    {
      Rprintf("i=%i\n",i);
    }
  
  
  /*---while文---*/
  i = 0;
  Rprintf("loop: increasing\n");
  while(i < 3)
    {
      Rprintf("i=%i\n",i);
      ++i;
    }
  
  i = 3;
  Rprintf("loop: decreasing\n");
  while(i > 0)
    {
      Rprintf("i=%i\n",i);
      --i;
    }
  
  
  
  /*---do-while文*/
  i = 0;
  Rprintf("loop: increasing\n");
  do
    {
      Rprintf("i=%i\n",i);
      ++i;
    }while(i < 3);
  
  i = 3;
  Rprintf("loop: decreasing\n");
  do
    {
      Rprintf("i=%i\n",i);
      --i;
    }while(i > 0);
}

/*** R
fSample5() 
*/
```

```
loop: increasing
i=0
i=1
i=2
loop: decreasing
i=3
i=2
i=1
loop: increasing
i=0
i=1
i=2
loop: decreasing
i=3
i=2
i=1
loop: increasing
i=0
i=1
i=2
loop: decreasing
i=3
i=2
i=1
```