{Rcpp}と{RcppGSL}を活用した数値積分の例

GSLはGNU Scientific Libraryの略で, 広い用途の数値計算向けのC言語のライブラリです。20年前から開発されており、まだリリースされ続けています。個人的な印象としては、めちゃくちゃチューニングされてはいないが、長年の開発のおかげで安定しており、マニュアル（htmlはこちら、pdfはこちら）も充実していて使いやすいイメージです。

この記事では{Rcpp}と{RcppGSL}パッケージを通してRからGSLを使う例（数値積分）を挙げます。このような{Rcpp}を使って数値計算を行う類似パッケージには{RcppNumerical}が挙げられます。

なお、RからC++を使うためのパッケージである{Rcpp}自体の解説はしません。{Rcpp}の使い方とメリットに関しては、最近非常によい本が出版されました。以下の本の24章を参照するとよいと思います。

パーフェクトR

インストール

GSLのインストール

Linuxの場合

僕の動作環境はRは3.3.1、gccは6.1.1、gslは2.3です。

パッケージマネージャのようなものからインストールすることもできますが、僕は諸事情でソースからインストールしました。Linuxで典型的な手順を踏めばインストールできます。例えば、公式ページからtar.gzをダウンロードして、解凍してconfigure、make、make installしてしばらく待ちます。

その後、.bashrcなどで「gslをインストールしたディレクトリ」を以下のように環境変数LD_LIBRARY_PATHに追加しておきます。

export LD_LIBRARY_PATH=".:「gslをインストールしたディレクトリ」/lib:$LD_LIBRARY_PATH"

Macの場合

持ってないので分かりません。すみません。

Windowsの場合

僕の動作環境はRは3.3.1、RtoolsはRtools33またはRtools34です。

最新版のGSLでなくてもよければ、Ripley先生が提供しているサイトから、Windows用にprebuildされたGSL（local323.zip）をダウンロードして、、、と思ったら、先週ぐらいから403 Forbiddenになっているようです。終了です。*1Ripley先生が退職されたあとも今まで残っていたのが奇跡なのかもしれません。Cygwinからだとインストールできるかもですが、Rtoolsとの連携もしんどそうで頑張る元気がありません。素直にLinuxかMacを使った方がよさそうです。

{Rcpp}と{RcppGSL}のインストール

いつも通りでOKです。すなわち、Rを起動してinstall.packages(c('Rcpp', 'RcppGSL'))でいけます。

テスト

テスト用のコードは作者のページにあるものを使います。C++ファイルとしてtest0.cppというファイル名で以下の内容で保存します。[[Rcpp::depends(RcppGSL)]]に注意です。

// [[Rcpp::depends(RcppGSL)]]
#include <RcppGSL.h>
#include <gsl/gsl_matrix.h>
#include <gsl/gsl_blas.h>

// [[Rcpp::export]]
Rcpp::NumericVector colNorm(const RcppGSL::Matrix & G) {
  int k = G.ncol();
  Rcpp::NumericVector n(k);
  for (int j = 0; j < k; j++) {
    RcppGSL::VectorView colview = gsl_matrix_const_column(G, j);
    n[j] = gsl_blas_dnrm2(colview);
  }
  return n;
}

さらに、別途Rファイルとしてtest0.Rというファイル名で以下の内容で保存します。

library(Rcpp)
library(RcppGSL)

sourceCpp('test0.cpp')

set.seed(123)
m <- matrix(rnorm(4), ncol=2)
print(colNorm(m)) #=> 0.6059 1.5603

• 4行目：ここで先ほど作成したC++ファイルへのパスを記入します。同じディレクトリに入っていれば上のままでOKです。

Rコンソールを起動してそのディレクトリに移動し、source('test0.R')を実行して無事に0.6059 1.5603が出力されれば完了です（小数点何位まで出力されるかはRのオプションの指定によって異なると思います）。

RcppGSL活用例：数値積分

Rに数値積分するintegrate関数があるのに、わざわざ自分で数値積分を書くバカがどこにおるの？と思われるかもしれませんが、近頃おじさんはRのintegrate関数が機能しない場面にしょっちゅう出くわします。例を挙げます。

ここで、とは値が与えられています。例えば、とし、をからの間のいくつか選んで算出してみます。数値積分の範囲はからまできちんと計算しなくても、からまでの積分でほぼ十分とします。

Rだけで算出する場合は以下になります。

sd <- 0.24
integrand <- function(x, mu) (0.2 + min(x, -0.3-x)) * dnorm(x, mean=mu, sd=sd)
g1 <- function(mu) integrate(integrand, mu-6*sd, mu+6*sd, mu)$value
g2 <- function(mu) pracma::romberg(integrand, a=mu-6*sd, b=mu+6*sd, mu=mu)$value

x <- seq(-3, 3, len=101)
m1 <- sapply(x, g1)
m2 <- sapply(x, g2)

g1はintegrate関数で、g2は比較ためのRomberg法で数値積分した結果です。{pracma}パッケージのromberg関数で実行できます。この結果を可視化すると以下になります。

値がだいぶずれていますが、romberg関数の方が正しい値です。それならばromberg関数をいつも使えばいいじゃんと思うかもしれませんが、計算に時間がかかるのと、ここでは挙げませんがromberg関数だとかえって数値計算がうまくいかない場合もあります。そこで、{RcppGSL}パッケージからで数値積分します。

テストの場合と同様に、C++ファイルとしてgsl_cquad.cppというファイル名で以下の内容で保存します。

// [[Rcpp::depends(RcppGSL)]]
#include <RcppGSL.h>
#include <gsl/gsl_integration.h>
#include <gsl/gsl_randist.h>

struct my_f_params { double mu; double sd; };

double integrand(double x, void *p) {
  struct my_f_params * params = (struct my_f_params *)p;
  double mu = params->mu;
  double sd = params->sd;
  double res = (0.2 + fmin(x, -0.3-x)) * gsl_ran_gaussian_pdf(x-mu, sd);
  return(res);
}

// [[Rcpp::export]]
Rcpp::NumericVector integrate_my_cquad(Rcpp::NumericVector xVec, double sd) {
  Rcpp::NumericVector resVec(xVec.size());
  gsl_integration_cquad_workspace *work = gsl_integration_cquad_workspace_alloc(1e2);

  double result(0.0);
  double abserr(0.0);
  size_t nevals(1e2);
  double mu(xVec[0]);

  gsl_function F;
  F.function = &integrand;
  struct my_f_params params = { mu, sd };

  for (int i=0; i < xVec.size(); ++i) {
    mu = xVec[i];
    params.mu = mu;
    F.params = &params;
    int res = gsl_integration_cquad(&F, mu-6*sd, mu+6*sd, 1e-8, 1e-6, work, &result, &abserr, &nevals);
    resVec[i] = result;
  }

  gsl_integration_cquad_workspace_free(work);
  return resVec;
}

• 3行目： GSLの数値積分の機能を使用する場合に必要となるヘッダーファイルです。
• 4行目： GSLの確率分布を使用する場合に必要となるヘッダーファイルです。
• 8～14行目： 被積分関数です。xが積分変数で、pはその他の変数（パラメータ）です。RcppGSLで数値積分する場合、被積分関数はR側で作って渡すのではなく、C++側で定義しないと高速にならないので注意です。なお、gsl_ran_gaussian_pdfは正規分布の確率密度関数です。
• 6行目：　いま説明したpに値を渡すための構造体です。
• 17～40行目：　数値積分を実行する自作関数を定義しています。RのvectorがxVecに渡され、そのxVecの各値ごとに数値積分を繰り返し実行して、その結果をRのvectorとなるresVecで返します。数値積分にもいろいろ手法があります。特異点がある場合にRのintegrate関数がこける印象があることから、QAGSとCQUADを検討した結果、CQUADの方が安定で速いケースが多かったので、この記事ではCQUADを使います。
• 34行目：　gsl_integration_cquadは数値積分をCQUADで実行する関数です。引数の説明についてはマニュアルを参照してください。引数の一つであるworkは作業領域になります。19行目のgsl_integration_cquad_workspace関数でその領域を確保し、38行目のgsl_integration_cquad_workspace_free関数でその領域を開放します。

C++ファイルで定義したintegrate_my_cquadを使うRスクリプトは以下になります。

library(Rcpp)
library(RcppGSL)

sourceCpp('gsl_cquad.cpp')

sd <- 0.24
x <- seq(-3, 3, len=101)
m3 <- integrate_my_cquad(x, sd)

計算の実行速度は非常に速く、また安定で、結果も上記のRomberg法とほぼ一致します。僕の実行環境で1000回ほどintegrate_my_cquad関数と同じ処理をして比較したところ、GSL版はRのromberg関数の約150倍、Rのintegrate関数の約1.5倍高速でした。

数値積分の簡単な使い方はこのサイトを参考にしました。

*1:portingのページはかろうじて残っていますが、、

MCMCサンプルを{dplyr}で操る

RからStanやJAGSを実行して得られるMCMCサンプルは、一般的に　iterationの数×chainの数×パラメータの次元　のようなオブジェクトとなっており、凝った操作をしようとするとかなりややこしいです。

『StanとRでベイズ統計モデリング (Wonderful R)』のなかでは、複雑なデータ加工部分は場合によりけりなので深入りしないで、GitHub上でソースコードを提供しています。そこでは、ユーザが新しく覚えることをなるべく少なくするため、Rの標準的な関数であるapply関数群を使っていろいろ算出しています。しかし、apply関数群は慣れていない人には習得しづらい欠点があります。

一方で、Rのデータ加工パッケージとして、%>%によるパイプ処理・{dplyr}パッケージ・{tidyr}パッケージがここ最近よく使われており、僕も重い腰を上げてやっと使い始めたのですが、これが凄く使いやすい。%>%、select、filter、mutate、group_by、summarize、*_join、pivot_longer、pivot_widerだけをまずは覚えればほとんど不自由しませんでした。これらがないともう他の言語に移れないレベルです。これらのパッケージの練習のおかげで、ややこしいMCMCサンプルの処理についても、こんな感じでやれば毎回ウンウン唸らずに統一的に操作できそうかなぁ、というところまで来ましたので簡単にメモします。

＊　＊　＊

手始めに以下の図を描いてみます。

この図はパラメータごとにMCMCサンプルの中央値と95%CIを表示した図です。{ggmcmc}パッケージや{bayesplot}パッケージに含まれる関数を使うと一撃で描くこともできます。しかし、練習のため自分で算出して作図します。

library(rstan)
library(ggmcmc)
library(dplyr)

data <- list(J=8, y=c(28,  8, -3,  7, -1,  1, 18, 12), sigma=c(15, 10, 16, 11,  9, 11, 10, 18))
model_code <- readr::read_file(url('https://raw.githubusercontent.com/wiki/stan-dev/rstan/8schools.stan'))
fit <- stan(model_code=model_code, data=data, seed=1234)

d_mcmc <- ggs(fit)
d_qua <- d_mcmc %>%
  filter(grepl('^theta\\[\\d+\\]$', Parameter)) %>%
  group_by(Parameter) %>%
  summarize(`2.5%` = quantile(value, probs=.025),
            `50%`  = quantile(value, probs=.5),
            `97.5%`= quantile(value, probs=.975)) %>% ungroup()

p <- ggplot() +
  geom_pointrange(data=d_qua, mapping=aes(x=forcats::fct_rev(Parameter), y=`50%`, ymin=`2.5%`, ymax=`97.5%`)) +
  coord_flip() +
  labs(x='Parameter', y='Value')
ggsave(p, file='fig1.png', dpi=300, w=4, h=3)

• 6行目：Web上から8schools.stanを読み込んで文字列としています。RStanの公式ページの例題で使われているモデルファイルです。
• 7行目：stan関数はmodel_code引数でモデルを書いた文字列も指定できます（本来はファイル名を直接指定できればよかったのですがよくわかりませんでした）。
• 9行目：{ggmcmc}パッケージのggs関数でtidyなデータにしておきます。tidyなデータについては西原さんの記事を参照。{ggmcmc}パッケージに含まれる関数を使うとd_mcmcからいろいろな図が描けます。詳しくは『StanとRでベイズ統計モデリング (Wonderful R)』の4章に書きましたので読んでいただけるとうれしいです。

なお、d_mcmcは以下のようなデータフレームになります。

> d_mcmc
# A tibble: 72,000 × 4
   Iteration Chain Parameter   value
       <dbl> <int>    <fctr>   <dbl>
1          1     1        mu -1.3476
2          2     1        mu -0.9601
3          3     1        mu  7.0919
4          4     1        mu 15.0782
5          5     1        mu 20.0110
6          6     1        mu 20.4483
7          7     1        mu 13.0249
8          8     1        mu 11.8232
9          9     1        mu 15.9213
10        10     1        mu 17.9294
# ... with 71,990 more rows

• 11行目：まずはモデルに含まれるtheta[数字]というパラメータだけ残しています。grepl関数でパラメータ名がマッチするか判定する際に正規表現を使う必要があります。ここが正規表現に慣れていない人は少し厳しいかもしれません。
• 12～15行目：{dplyr}パッケージの典型的な使い方です。Parameter列ごとに要約量を算出します。列名が数字で始まる場合はバッククォートで囲む必要があります。1つずつ分位点を算出するのではなく、do関数で一行で算出する方法もあるのですが、分かりにくく、現状issueとして検討中のようです（ここ→ここ）。
• 18行目：ggplot2でcoord_flipすると下から上に向かってfactorが並びますので、{forcats}パッケージのfct_rev関数で逆順にしています。

＊　＊　＊

次に以下の図を描いてみます。久保本11章の図に相当します。

library(rstan)
library(ggmcmc)
library(dplyr)

Y <- read.csv(url('https://raw.githubusercontent.com/MatsuuraKentaro/RStanBook/master/chap12/input/data-kubo11a.txt'))$Y
I <- length(Y)
d <- data.frame(X=1:I, Y=Y)
data <- list(I=I, Y=Y)
model_code <- readr::read_file(url('https://raw.githubusercontent.com/MatsuuraKentaro/RStanBook/master/chap12/model/model12-11.stan'))
fit <- stan(model_code=model_code, data=data, seed=1234)

d_mcmc <- ggs(fit)
d_qua <- d_mcmc %>%
  filter(grepl('^Y_mean\\[\\d+\\]$', Parameter)) %>%
  tidyr::separate(Parameter, into=c('Parameter', 'x'), sep='[\\[\\]]', convert=TRUE) %>%
  group_by(Parameter, x) %>%
  summarize(`2.5%` = quantile(value, probs=.025),
            `10%`  = quantile(value, probs=.1),
            `50%`  = quantile(value, probs=.5),
            `90%`  = quantile(value, probs=.9),
            `97.5%`= quantile(value, probs=.975)) %>% ungroup()

p <- ggplot() +
  geom_ribbon(data=d_qua, mapping=aes(x=x, ymin=`2.5%`, ymax=`97.5%`), alpha=1/6) +
  geom_ribbon(data=d_qua, mapping=aes(x=x, ymin=`10%`,  ymax=`90%`),   alpha=2/6) +
  geom_line(data=d_qua, mapping=aes(x=x, y=`50%`)) +
  geom_point(data=d, aes(x=X, y=Y), shape=1, size=2) +
  labs(x='i', y='Y[i]') +
  ylim(0, 22)
ggsave(p, file='fig2.png', dpi=300, w=4, h=3)

• 15行目：{tidyr}パッケージのseparate関数を使って、Y_mean[20]をY_meanと20という2つの列に分解しています。
• 16行目：あとは集計の単位であるgroup_byの単位が場面によって多少変わるぐらいで、特に悩まずに色々な量が算出できます。

この記事ではMCMCのChainごとに何かを算出することは取り上げませんでしたが、ggs関数で作ったd_mcmcはChain列も含んでいますので自由自在です。

2017.07.16 追記

vector[D] mu[T]のようにD次元vectorがT（時点の数）個並べたような配列に対し、各t,dのmuの値の分位点を算出したい場合は以下のようにしました。tidyr::separateでまず角括弧の切れ目で分けておいてから、カンマで切るのがポイントです。

d_qua <- d_mcmc %>%
  filter(grepl('^mu\\[\\d+,\\d+\\]$', Parameter)) %>%
  tidyr::separate(Parameter, into=c('Parameter', 't'), sep='[\\[\\]]') %>%
  tidyr::separate(t, into=c('t', 'd'), sep=',', convert=TRUE) %>%
  group_by(Parameter, t, d) %>%
  summarize(`2.5%` = quantile(value, probs=.025),
            `10%`  = quantile(value, probs=.1),
            `50%`  = quantile(value, probs=.5),
            `90%`  = quantile(value, probs=.9),
            `97.5%`= quantile(value, probs=.975)) %>% ungroup()

2次元以上のスポット検出を行う統計モデル

1次元の場合の変化点検出は以下の記事で扱いました。

• 状態空間モデルでシステムノイズに非ガウス分布（1次元の変化点検出） - StatModeling Memorandum
• 二つの時系列データの間に「差」があるか判断するには - StatModeling Memorandum

変化点検出のポイントは状態空間モデルのシステムモデルに「ほとんどの場合は0の近くの値を生成するが，まれにとても大きな値を生成する」という性質を持つ分布を使うことでした。そこで、上記の例ではコーシー分布を使いました。しかし、コーシー分布はかなり厄介な分布なので逆関数法を利用したコーディングをしないとうまくサンプリングできません。そのため2次元以上の場合の同様のモデル（マルコフ場モデル）に拡張することができません。2次元の場合のマルコフ場モデルについては以下の記事で簡単に扱いました（少しコードが古いです）。詳しくは書籍「StanとRでベイズ統計モデリング」の12章を参照してください。

• 久保緑本11章のマルコフ場モデル（空間構造のあるベイズモデル） - StatModeling Memorandum

それでも2次元以上のスポット検出を行いたくて色々なモデルを試しました。コーシー分布＋実装の工夫、線過程（ライン過程）、混合正規分布、ベルヌーイ分布と正規分布の混合分布、両側指数分布（ラプラス分布）を使ったモデルなどなど。しかし、少なくとも私の環境においてはこれらのモデルだと、推定がうまくいかない場合や、推定できてもスポットのようにシャープではなくボワッと境界があいまいに推定されてしまう場合が多かったです。その後も2年にわたり粘り続けて合計150個弱のモデルを試した結果、一つの解決策を見つけましたのでこの記事を書きます。

UBN分布

アイデアは「ほとんどの場合は0の近くの値を生成するが，まれにとても大きな値を生成する」という性質を工夫して表現することです。そこで、状態空間モデルのシステムモデルに相当する部分に以下の分布を使いました。

この分布は]の範囲の値を生成する、切断正規分布に一様分布のゲタをはかせたような分布です。ここで、は一様分布の範囲を決めるパラメータ、は一様分布のゲタの高さを決めるパラメータ、は正規分布の平均と標準偏差、は正規化項です。この分布は一様分布と切断正規分布の混合分布なのですが、ベルヌーイ分布とポアソン分布の混合分布をZero-inflated Poisson分布と分かりやすく呼ぶのになぞらえて、この記事では上記の分布をUniformly Boosted Normal Distributionと呼び、以下では略してUBN分布と呼びます。思いつけばコロンブスの卵ですが、UBN分布はコーシー分布の利点である裾の長さを保ったまま、サンプリングの非効率さを解消します。例えば、,,,,の場合にははほぼ2となり、以下のような確率密度になります。

なお、Stanの計算に必要な対数尤度は以下のように式変形して表現することで高速化できます。

ここでlog1p_exp関数とnormal_lpdf関数は計算の高速化するためのStanの便利関数です。

2次元のスポット検出例 その1

30x30の大きさで中央部分にシグナルがあり、ノイズを加えたデータYを考えます。図にすると以下です。

上下左右のマスとつながるマルコフ場モデルを、UBN分布を利用して実装したStanコードは以下になります。

• 5行目： 上記で説明したUBN分布のを固定値Uで与えています。
• 9行目：　muにノイズがのってデータYになります。設定した下限と上限がUBN分布のとにそれぞれ対応しています。
• 17,20行目： 17行目がマルコフ場モデルの左右のつながり、20行目が上下のつながりに対応します。ここではUを固定値として与えていること、またs_muの値が高々1.5ぐらいと考えられることから、UBN分布のはほぼ一定の定数です。Stanでは対数確率の偏微分が重要であり、微分で消える定数項を無視して構いません。そのためを無視してコーディングしています。
• 21行目： ノイズをのせている部分です。状態空間モデルの観測モデルに相当します。
• 22,23行目： なくても構いませんが、変分ベイズが安定になることがあるのでs_muとs_Yに弱情報事前分布を設定しています。

このモデルは局所最適解が多くMCMCの一種であるNUTSではうまく推定できませんが、自動変分ベイズだと高速に推定できるモデルとなっています。データを生成し、推定を実行するRコードは以下になります。

• 4～9行目：データ生成部分です。
• 12行目：ここではUは0.1としました。少し値を変えて実行してみて結果が頑健かチェックするとよいでしょう。

変分ベイズが高速なため、計算時間はSurface Pro 3（core i5）でも10秒以下です。推定結果は以下です。以下の図はmuの乱数サンプルの中央値を2次元でプロットした図です。

以下の図はRowの真ん中（ここでは15）に固定してスライスし、横軸をColumn・縦軸をmuにしてプロットした図です。薄い灰帯は乱数サンプルから計算した95%ベイズ信頼区間、濃い灰帯は80%ベイズ信頼区間、青線は中央値、赤線はデータです。

シグナルのスポットを検出できているのがわかるかと思います。上記のRコードと結果はシグナルの大きさがノイズのSDの3倍となっている場合ですが、シグナルの大きさがSDの2.5倍となっている場合の結果は以下になります。徐々にカドから崩れていきます。

このシグナルが正方形の場合、少なくともシグナルがノイズのSDの2.3倍ぐらいないとスポットとして綺麗に検出されないようでした。

2次元のスポット検出例 その2

画素数が300x600の場合も試してみました。データ例を図にすると以下です。

モデルは同じままで、実行するRコードもほぼ同じです（データ生成部分だけわずかに異なります）。推定に要する時間は5分ぐらいでした。シグナルがノイズのSDの3倍の場合の推定結果は以下になりました。図の凡例は上記の場合と同じです。丸いスポットの検出は比較的得意ですが、尖っている箇所はやはり難しいようです。

この場合、少なくともシグナルがノイズのSDの3倍ぐらいないとスポットとして綺麗に検出されないようでした。

試した範囲では、この方法によって地図上でも3次元でもスポット検出できそうです。Enjoy!