Stanの関数を使ってRを拡張して高速化する

（C++に自動で変換される）Stanの関数を使ってRを拡張できる機能が、Stan/RStanの2.16で実装開始されて2.17でほぼ完成しました。Rを高速化するためにC++（とRcpp）はあまり書きたくないけれど、Stanの関数なら書いてもいいよという僕得な機能です。この記事ではその方法を簡単に紹介します。

元にした資料はRStanの開発者であるBenさんがStanCon2018で発表したこちらの資料です。

ここでは例として、以下の2つの関数をRで使えるようにしましょう。

• 1) 機械学習分野でおなじみのlog_sum_exp関数
  • 引数はN個の正の実数
• 2) データにemax modelという曲線をあてはめた場合の対数尤度を返す関数
  • 引数はデータ（N個のXとYのペア）とパラメータの値

手順は簡単で以下だけです。

• functionsブロックだけ書いたstanファイルを用意する
• R側でrstan::expose_stan_functions()を呼び出す

具体的には、以下の内容を例えばmy_functions.stanという名前で保存します。

functions {
  real stan_log_sum_exp(vector x) {
    return(log_sum_exp(x));
  }

  real emax_model_ll(vector Y, vector X, real e0, real emax, real ed50, real sigma) {
    int N = num_elements(Y);
    vector[N] mu = e0 + emax*X ./ (ed50 + X);
    real ll = normal_lpdf(Y | mu, sigma);
    return(ll);
  }
}

• 3行目：log_sum_exp関数の方はすでにStanに用意されているものを使います。
• 6行目：接尾の_llはlog likelihoodの意味です。
• 7行目：たしかStan 2.14ぐらいから変数の宣言と定義を同時にできるようになりました。ここではデータポイントの数NをYから動的に取得しています。

次にR側でrstan::expose_stan_functions()を呼び出します。

rstan::expose_stan_functions('my_functions.stan')

するとfunctionsブロックに入っていた関数は全てR側で使えるようになります。

stan_log_sum_exp(c(0.1, 0.05, 0.8))
emax_model_ll(Y=c(1,10,12,24), X=c(0,5,10,10), e0=0.8, emax=30, ed50=5.1, sigma=3)

StanのコードはC++のコードに変換されるのでfor文はRと比べて非常に高速です。また、Stanのマニュアル（特にBuilt-In Functionsの章以降）を見ていただくとよいのですが、Stanはマニアックな分布の対数尤度を求める関数や、行列演算の関数（Eigen由来）や、log1m_expなどの指数対数を安定に計算する関数が豊富です。

これらのStanの長所を非常に簡単にRのコーディングに取り込めるのは最高です。まあ、少しばかりStanの関数やデータ型に慣れている必要がありますが。

Enjoy!