本記事では、BayesRTMB が内部で利用している
RTMB
パッケージの仕組みと、モデルの推定がどのようなアルゴリズムで行われているかをモデラー向けに解説します。
1. RTMBとは何か?
RTMB は、C++で書かれた強力な自動微分ライブラリ
TMB (Template Model Builder)
を、Rの構文から直接利用できるようにしたパッケージです。
テーピング(Taping)の仕組み
RTMB
の最大の特徴は「テーピング」と呼ばれるプロセスにあります。rtmb_model()
を実行した際、Rのコードは一度だけ実行され、その計算過程が「計算グラフ(テープ)」として記録されます。
-
実行: Rの関数(
rtmb_codeで定義したもの)が、「AD(自動微分)タイプ」のオブジェクトを引数として実行されます。 -
記録: 四則演算や数学関数(
log,exp,sinなど)が呼び出されるたびに、その操作がテープに記録されます。 - コンパイル: 記録されたテープは、裏側のC++エンジンで高速に評価可能な形式に変換されます。
この仕組みにより、ユーザーはRの書きやすさを維持したまま、C++並みの速度で勾配(導関数)やヘッセ行列を計算できるようになります。Stanのようにコードを修正するたびにコンパイルをし直す必要がなく、すぐにMCMCやMAP推定を実行できます。
2. 推定のアルゴリズム
BayesRTMB の optimize() (MAP推定) や
sample() (MCMC)
は、Stanと同様の高度なアルゴリズムを採用しています。
パラメータ空間の変換
多くの統計モデルでは、パラメータに制約(例:標準偏差は正、確率は0から1の間、相関行列の性質など)があります。しかし、勾配ベースの最適化アルゴリズムや HMC (Hamiltonian Monte Carlo) は、「制約のない実数空間 ()」で最も効率的に動作します。
BayesRTMB は、以下のプロセスを自動的に行います。
-
制約付き空間
():
ユーザーが定義したパラメータ(例:
lower = 0のsigma)。 -
変換
():
制約のない実数空間
()
への写像。
- 正数制約:
- 上限下限制約:
- 実数空間 () での推定: 最適化やサンプリングはこの空間で行われます。
ヤコビアン(Jacobian)の調整
変数変換を行った場合、変換後の空間での確率密度を正しく保つために、変換のヤコビアンの行列式の絶対値の対数を加算する必要があります。
BayesRTMB は、rtmb_model()
を作成する際にこのヤコビアン調整を自動的にモデルの対数密度に加算します。これにより、ユーザーはヤコビアンを意識することなく、Stanと同じように正しい事後分布からのサンプリングが可能になります。
3. RTMBコードを書く際の注意点
RTMBはRのように見えますが、内部は「計算グラフの記録」を行っているため、通常のRとは異なる制限があります。
条件分岐の制限
計算グラフは一度しか記録されないため、パラメータの値によって計算パスが変わる
if 文は、期待通りに動作しません。
-
誤り:
if (sigma > 1) { ... } else { ... } -
正しい:
RTMB::adm_if,RTMB::adm_ifelseを使うか、滑らかな関数(abs,plogisなど)で代用する。