書評.この広告は、90日以上更新していないブログに表示しています。広い用途の数値計算向けのC言語のライブラリです。20年前から開発されており、まだリリースされ続けています。個人的な印象としては、めちゃくちゃチューニングされてはいないが、長年の開発のおかげで安定しており、マニュアル(htmlはこちら、pdfはこちら)も充実していて使いやすいイメージです。この記事では{Rcpp}と{RcppGSL}パッケージを通してRからGSLを使う例(数値積分)を挙げます。このような{Rcpp}を使って数値計算を行う類似パッケージには{RcppNumerical}が挙げられます。RからC++を使うためのパッケージである{Rcpp}自体の解説はしません。{Rcpp}の使い方とメリットに関しては、最近非常によい本が出版されました。以下の本の24章を参照するとよいと思います。僕の動作環境はRは3.3.1、gccは6.1.1、gslは2.3です。パッケージマネージャのようなものからインストールすることもできますが、僕は諸事情でソースからインストールしました。Linuxで典型的な手順を踏めばインストールできます。installしてしばらく待ちます。その後、.bashrcなどで「gslをインストールしたディレクトリ」を以下のように環境変数LD_LIBRARY_PATHに追加しておきます。持ってないので分かりません。僕の動作環境はRは3.3.1、RtoolsはRtools33またはRtools34です。Forbiddenになっているようです。*1Ripley先生が退職されたあとも今まで残っていたのが奇跡なのかもしれません。Cygwinからだとインストールできるかもですが、Rtoolsとの連携もしんどそうで頑張る元気がありません。素直にLinuxかMacを使った方がよさそうです。{Rcpp}と{RcppGSL}のインストールいつも通りでOKです。'RcppGSL'))でいけます。テストテスト用のコードは作者のページにあるものを使います。C++ファイルとしてtest0.cppというファイル名で以下の内容で保存します。[[Rcpp::depends(RcppGSL)]]に注意です。n;}さらに、別途Rファイルとしてtest0.Rというファイル名で以下の内容で保存します。1.56034行目:ここで先ほど作成したC++ファイルへのパスを記入します。同じディレクトリに入っていれば上のままでOKです。1.5603が出力されれば完了です(小数点何位まで出力されるかはRのオプションの指定によって異なると思います)。RcppGSL活用例:数値積分Rに数値積分するintegrate関数があるのに、わざわざ自分で数値積分を書くバカがどこにおるの?と思われるかもしれませんが、近頃おじさんはRのintegrate関数が機能しない場面にしょっちゅう出くわします。値が与えられています。間のいくつか選んで算出してみます。積分でほぼ十分とします。Rだけで算出する場合は以下になります。g2)g1はintegrate関数で、g2は比較ためのRomberg法で数値積分した結果です。{pracma}パッケージのromberg関数で実行できます。この結果を可視化すると以下になります。値がだいぶずれていますが、romberg関数の方が正しい値です。それならばromberg関数をいつも使えばいいじゃんと思うかもしれませんが、計算に時間がかかるのと、ここでは挙げませんがromberg関数だとかえって数値計算がうまくいかない場合もあります。そこで、{RcppGSL}パッケージからで数値積分します。テストの場合と同様に、C++ファイルとしてgsl_cquad.cppというファイル名で以下の内容で保存します。GSLの数値積分の機能を使用する場合に必要となるヘッダーファイルです。GSLの確率分布を使用する場合に必要となるヘッダーファイルです。xが積分変数で、pはその他の変数(パラメータ)です。RcppGSLで数値積分する場合、被積分関数はR側で作って渡すのではなく、C++側で定義しないと高速にならないので注意です。gsl_ran_gaussian_pdfは正規分布の確率密度関数です。いま説明したpに値を渡すための構造体です。数値積分を実行する自作関数を定義しています。RのvectorがxVecに渡され、そのxVecの各値ごとに数値積分を繰り返し実行して、その結果をRのvectorとなるresVecで返します。数値積分にもいろいろ手法があります。特異点がある場合にRのintegrate関数がこける印象があることから、QAGSとCQUADを検討した結果、CQUADの方が安定で速いケースが多かったので、この記事ではCQUADを使います。gsl_integration_cquadは数値積分をCQUADで実行する関数です。引数の説明についてはマニュアルを参照してください。引数の一つであるworkは作業領域になります。19行目のgsl_integration_cquad_workspace関数でその領域を確保し、38行目のgsl_integration_cquad_workspace_free関数でその領域を開放します。C++ファイルで定義したintegrate_my_cquadを使うRスクリプトは以下になります。sd)計算の実行速度は非常に速く、また安定で、結果も上記のRomberg法とほぼ一致します。僕の実行環境で1000回ほどintegrate_my_cquad関数と同じ処理をして比較したところ、GSL版はRのromberg関数の約150倍、Rのintegrate関数の約1.5倍高速でした。数値積分の簡単な使い方はこのサイトを参考にしました。高速化&フルベイズ前回の記事のスピードアップをします。(2)月別アーカイブはてなブログをはじめよう!StatModelingさんは、はてなブログを使っています。ブログを報告するスターをつけました引用をストックしましたストック一覧を見る閉じる引用するにはまずログインしてくださいログイン閉じる引用をストックできませんでした。再度お試しください閉じる限定公開記事のため引用できません。