目的別クイックナビ

1.1 使い分けの目安

迷った場合は、まず次の表を見ます。

2.1 モデル定義

2.2 ラッパー関数

2.3 推定メソッド

sample() が標準的なベイズ推論です。optimize() と variational() は速い確認や初期値探索に使えますが、不確実性評価やBayes factorではMCMCを優先します。

3.1 estimate(), EAP(), MAP()

estimate() は、パラメータ、変換量、生成量を取り出すための共通APIです。

# primary parametersだけを取り出す
fit$estimate(pars = "parameters")

# transformブロックで作った量を取り出す
fit$estimate(pars = "transform")

# generateブロックで作った量を取り出す
fit$estimate(pars = "generate")

# すべて取り出す
fit$estimate(pars = "all")

pars = "parameters" は、parameters ブロックで宣言されたprimary parameterだけを返します。推定結果を初期値や固定値として再利用したいときは、通常これを使います。

est <- fit$estimate(pars = "parameters", drop = FALSE)

drop = TRUE がデフォルトです。選択されたパラメータが1つだけの場合は、リストではなくベクトル、行列、配列そのものを返します。常にリストで受け取りたい場合は drop = FALSE を指定します。

# bだけならベクトルで返る
b_eap <- fit$EAP("b")

# bだけでもlistとして返る
b_eap_list <- fit$EAP("b", drop = FALSE)

EAP() は estimate(type = "EAP") のショートカットです。MAP() は周辺MAPを返すのがデフォルトです。

fit$EAP(pars = "parameters")
fit$MAP(pars = "parameters")
fit$MAP(pars = "parameters", type = "joint")

3.2 parsとcomponent

estimate() では、pars と component を使って対象を絞れます。

# bとsigmaだけ
fit$estimate(pars = c("b", "sigma"))

# transformだけ
fit$estimate(component = "transform")

# generateだけ
fit$estimate(component = "generate")

# theta以外
fit$estimate(pars = "-theta")

MCMC/VBでは、pars = NULL の場合、デフォルトではprimary parametersが返ります。Classic/MAPでは、推定済みの固定効果、変換量、生成量を含む結果が返ります。

4.1 主要メソッド

4.2 draws()

# 全ての固定パラメータ、transform、generateを取得
draws_all <- fit_mcmc$draws()

# bだけ取得
draws_b <- fit_mcmc$draws("b")

# random effectsも含める
draws_with_random <- fit_mcmc$draws(inc_random = TRUE)

# transformやgenerateを含めない
draws_par <- fit_mcmc$draws(
  inc_transform = FALSE,
  inc_generate = FALSE
)

draws() の返り値は3次元配列です。第3次元の名前がパラメータ名です。

dim(draws_b)
dimnames(draws_b)[[3]]

4.3 summary()とdiagnose()

fit_mcmc$summary()
fit_mcmc$summary("b")
fit_mcmc$diagnose()

summary() は、mean、sd、map、q2.5、q97.5、ess_bulk、ess_tail、rhat を表示します。

diagnose() では、次のような情報を確認します。

• retained draws
• finite lp
• mean acceptance
• maximum treedepth
• divergent transitions
• mean/max leapfrog steps
• metric type
• warmup summary
• dense metric condition number
• R-hat
• ESS
• pointwise log_lik の有無

MCMCの最終確認では、R-hat、ESS、divergence、treedepth、acceptanceをセットで確認します。

d <- fit_mcmc$diagnose()
print(d)

よくある警告と、最初に試すことは次の通りです。

4.4 transformed_draws()とgenerated_quantities()

transform や generate の計算を、推定後にposterior drawごとに追加できます。

fit_mcmc$transformed_draws()

fit_mcmc$generated_quantities({
  y_rep <- rnorm(length(Y), mu, sigma)
})

generated_quantities() には、rtmb_code(generate = { ... }) または { ... } を渡せます。

gq_code <- rtmb_code(
  generate = {
    log_lik <- normal_lpdf(Y, mu, sigma, sum = FALSE)
  }
)

fit_mcmc$generated_quantities(gq_code)
fit_mcmc$WAIC()

progress = "message" を指定すると、Jobや並列実行でも進捗メッセージとして表示されます。

fit_mcmc$generated_quantities(gq_code, progress = "message")

4.5 bridge samplingとBayes factor

MCMC fitでは、bridge samplingで周辺尤度を推定できます。

logml <- fit_mcmc$bridgesampling(
  method = "normal",
  use_neff = TRUE
)

logml
attr(logml, "error")
attr(logml, "ess")

attr(logml, "error") はbridge samplingの推定誤差の目安、attr(logml, "ess") はbridge weightに基づく有効サンプルサイズの目安です。誤差が大きい、またはESSが小さい場合は、sampling や chains を増やす、method を見直す、posterior drawの混合を確認する、という順に疑います。

Bayes factorは、周辺尤度の比として計算されます。

bf <- fit_mcmc$bayes_factor(
  fixed = list("b[talk]" = 0),
  sampling = 4000,
  chains = 4
)

print(bf)

fixed = list(...) を指定すると、指定したパラメータを固定した比較モデルを内部で作り、そのモデルもMCMCで推定してからBayes factorを計算します。Bayes factorの解釈と診断は8.1節にもまとめています。

すでに比較モデルを推定している場合は、comparison_fit を渡せます。

mdl_null <- fit_mcmc$model$fixed_model(
  fixed = list("b[talk]" = 0)
)
fit_null <- mdl_null$sample(chains = 4, sampling = 4000)

bf <- fit_mcmc$bayes_factor(
  comparison_fit = fit_null
)

Bayes factorはpriorに敏感です。モデル比較を目的にする場合は、prior_flat() ではなく、分析目的に合ったproper priorを置きます。

4.6 WAIC

WAICには、観測ごとのlog likelihoodが必要です。自作モデルでは generate ブロックに log_lik を保存します。モデル比較としてのWAICの使い方は8.2節にまとめています。

df <- debate

code <- rtmb_code(
  setup = {
    Y <- sat
  },
  parameters = {
    mu <- Dim()
    sigma <- Dim(lower = 0)
  },
  model = {
    Y ~ normal(mu, sigma)
  },
  generate = {
    log_lik <- normal_lpdf(Y, mu, sigma, sum = FALSE)
  }
)

fit <- rtmb_model(data = df, code = code)$sample()
fit$WAIC()

ラッパー関数では、対応している場合に WAIC = TRUE を指定します。

mdl <- rtmb_lm(
  sat ~ talk * perf,
  data = debate,
  prior = prior_normal(),
  WAIC = TRUE
)

fit <- mdl$sample()
fit$WAIC()

5.1 主要メソッド

5.2 num_estimateとbest run

optimize(num_estimate = K) は、複数の初期値から最適化を行い、not convergedではないrunの中から最良の目的関数値を持つrunを選びます。最終的な fit_map$estimate() や fit_map$summary() は、選ばれたbest runに基づきます。

各runの状態は opt_history で確認します。

fit_map <- mdl$optimize(num_estimate = 8)
fit_map$opt_history
fit_map$diagnose()

局所解があり得るモデルでは、num_estimate を増やすことが有効です。

fit_mix <- mdl_mix$optimize(num_estimate = 20)
fit_mix$opt_history

5.3 optimize結果を初期値や固定値に使う

MAP推定結果は、MCMCの初期値としてよく使います。

fit_map <- mdl$optimize(num_estimate = 4)

fit_mcmc <- mdl$sample(
  init = fit_map$estimate(pars = "parameters", drop = FALSE)
)

一部のパラメータをMAP推定値で固定した制約付きモデルを作ることもできます。

est <- fit_map$estimate(pars = "parameters", drop = FALSE)

mdl_fixed <- mdl$fixed_model(
  fixed = list(
    sigma = est$sigma
  )
)

fit_fixed <- mdl_fixed$optimize()

ベクトルや行列パラメータの一部だけを固定する場合は、展開名を使います。

mdl_null <- mdl$fixed_model(
  fixed = list("b[talk]" = 0)
)

fit_null <- mdl_null$optimize()

5.4 profile()

profile() は、特定のパラメータについてprofile likelihood区間を計算します。

prof <- fit_map$profile(
  pars = c("b[talk]", "sigma"),
  level = 0.95
)

Laplace近似を使っているfitでは、内部的にrandom側のヤコビアン補正を含めてprofile用の目的関数を再構築します。

5.5 MAP_FitのWAIC

MAP_Fit$WAIC() は、se_method = "sampling" による不確実性伝播があり、かつ log_lik が保存されている場合に使えます。

mdl <- rtmb_lm(
  sat ~ talk * perf,
  data = debate,
  prior = prior_normal(),
  WAIC = TRUE
)

fit_map <- mdl$optimize(se_method = "sampling")
fit_map$WAIC()

MAPのWAICは近似的です。最終的な予測評価ではMCMCのWAICを優先します。

6.1 主要メソッド

6.2 num_estimateとbest estimate

variational(num_estimate = K) は、複数のVB推定を行い、ELBOが最も高いrunをbest estimateとして選びます。

VB_Fit では、summary()、estimate()、EAP()、MAP() は、デフォルトでbest estimateだけを使います。全estimateのdrawを確認したい場合は draws() を使います。

fit_vb$best_chain
fit_vb$ELBO
fit_vb$rel_obj_vals

# デフォルトではbest estimateのEAP
fit_vb$EAP(pars = "parameters")

# 特定のestimateを指定
fit_vb$EAP(pars = "parameters", chains = 2)

# ELBO上位2つを使う
fit_vb$EAP(pars = "parameters", best_chains = 2)

6.3 plot_elbo()による収束確認

VBでは、ELBOが上昇し、その後ほぼ水平になることが収束の目安です。

fit_vb$plot_elbo()
fit_vb$plot_elbo(tail_n = 1000)
fit_vb$plot_elbo(ests = "best", tail_n = 1000)

後半のELBOが明らかに上昇し続けている場合は、iter を増やします。複数のestimateでELBOが大きくばらつく場合は、num_estimate を増やすか、MCMCで確認します。

fit_vb$diagnose()

VBの収束確認では、次の順に見ます。

1. plot_elbo() の後半が水平に近いか。
2. rel_obj_vals が小さく、最適化が止まる理由として不自然でないか。
3. num_estimate 間でELBOが極端に違わないか。
4. diagnose() がVB objectiveやWAICについて警告していないか。

ELBOが水平であることは必要な目安ですが、十分条件ではありません。VBは近似なので、posterior uncertaintyが重要な結論ではMCMCで確認します。

6.4 VBのWAIC

VBでも log_lik があればWAICを計算できます。

mdl <- rtmb_lm(
  sat ~ talk * perf,
  data = debate,
  prior = prior_normal(),
  WAIC = TRUE
)

fit_vb <- mdl$variational()
fit_vb$WAIC()

ただし、VBのposteriorは近似であり、事後分散やWAICが楽観的になることがあります。モデル選択の判断では、MCMCでも確認します。

7.1 主要メソッド

Classic_Fit では EAP() と MAP() は使いません。点推定は estimate() を使います。

fit_classic$estimate()
fit_classic$summary()
fit_classic$diagnose()

Classic APIは、wrapper関数から作られる頻度主義的な出力を扱うためのものです。anova() はlm/glm/mixed modelや分割表など、モデルごとに定義された範囲で使えます。rtmb_ttest()、rtmb_corr()、rtmb_fa()、rtmb_irt() などでは、すべてのclassicメソッドが同じ意味で使えるわけではありません。lsmeans() は主にカテゴリカル因子を含む回帰系モデルで使う機能で、bootstrap() は現状では主にmediationでの利用を想定しています。

7.2 ANOVA

anova() は、モデルに応じてWald F test、Wald chi-square test、分割表検定、モデル比較を行います。

fit_classic$anova()
anova(fit_classic)

複数の Classic_Fit を渡すと、モデル比較になります。

fit0 <- rtmb_lm(sat ~ 1, data = debate, prior = prior_flat())$classic()
fit1 <- rtmb_lm(sat ~ talk + perf, data = debate, prior = prior_flat())$classic()

anova(fit0, fit1)

7.3 AIC, BIC, logLik

fit_classic$logLik()
fit_classic$AIC()
fit_classic$BIC()

AIC(fit0, fit1)
BIC(fit0, fit1)

AIC() と BIC() は、頻度主義的な情報量基準として使います。ベイズ的なモデル比較とは目的が異なります。

7.4 robust_se()

robust_se() は、lm/glm/mixed modelでrobust SEを計算します。

# 元のfitを変更せず、robust SE版のコピーを返す
fit_robust <- fit_classic$robust_se(type = "HC3")

# cluster-robust
fit_cluster <- fit_classic$robust_se(
  cluster = "group_id",
  type = "CR1"
)

# fit自体を更新
fit_classic$robust_se(type = "HC3", inplace = TRUE)

7.5 lsmeans()

lsmeans() は、カテゴリカル因子のmarginal meansやpairwise comparisonを計算します。

fit_classic$lsmeans("group")
fit_classic$lsmeans("group", pairwise = TRUE)
fit_classic$lsmeans(c("group", "condition"), pairwise = TRUE)

simple_effects() や conditional_effects() も、対応するfitで標準誤差や不確実性が利用できる場合に使えます。

conditional_effects(fit_classic, effect = "talk:perf")
simple_effects(fit_classic, effect = "talk:perf")

8.1 Bayes factor

Bayes factorは、2つのモデルの周辺尤度の比です。BayesRTMBでは、MCMC fitのbridge samplingに基づいて計算します。

fit_full <- mdl_full$sample(chains = 4, sampling = 4000)

bf <- fit_full$bayes_factor(
  fixed = list("b[talk]" = 0),
  chains = 4,
  sampling = 4000
)

fixed は、nested null modelを作るための指定です。fixed = list("b[talk]" = 0) は、b[talk] を0に固定したモデルと比較します。

Bayes factorを見るときは、値そのものだけでなく、bridge samplingの診断も確認します。

logml <- fit_full$bridgesampling()

attr(logml, "error")
attr(logml, "ess")

attr(logml, "error") が大きい場合は、周辺尤度の数値誤差がBayes factorの判断に影響する可能性があります。attr(logml, "ess") が小さい場合は、bridge weightの一部に依存した不安定な推定になっている可能性があります。bayes_factor() で error_threshold を指定すると、誤差が大きい場合に警告・停止しやすくできます。

bf <- fit_full$bayes_factor(
  fixed = list("b[talk]" = 0),
  chains = 4,
  sampling = 8000,
  error_threshold = 0.05
)

Bayes factorはpriorに敏感です。比較したい仮説に対応するproper priorを事前に決めておきます。prior_flat() のような不適切事前分布は、周辺尤度やBayes factorの比較には向きません。

8.2 WAIC

WAICは、事後分布に基づく予測評価です。観測ごとの log_lik が必要です。

mdl <- rtmb_lm(
  sat ~ talk * perf,
  data = debate,
  prior = prior_normal(),
  WAIC = TRUE
)

fit1 <- mdl$sample()
fit1$WAIC()

自作モデルでは generate ブロックに log_lik を置きます。

generate = {
  log_lik <- normal_lpdf(Y, mu, sigma, sum = FALSE)
}

sum = FALSE にして、観測ごとのlog likelihoodを返します。

8.3 AIC/BIC

AIC/BICは Classic_Fit に対して使います。

fit0 <- mdl0$classic()
fit1 <- mdl1$classic()

AIC(fit0, fit1)
BIC(fit0, fit1)
anova(fit0, fit1)

8.4 ELBO

VBではELBOが最適化対象です。ELBOが大きいほど、近似事後分布の目的関数上はよい解です。

fit_vb$ELBO
fit_vb$plot_elbo(tail_n = 1000)

ELBOは異なる推定法間の一般的なモデル比較指標ではありません。VBの収束確認や、同じモデル・同じ設定内のrun選択に使います。

9.1 固定するパラメータ名の確認

fixed には、モデル内部で使われるパラメータ名を渡します。ベクトルや行列の一部を固定する場合は、"b[talk]" のような展開名を使います。

名前を確認するには、次の方法が便利です。

# summaryに表示される名前を見る
fit_mcmc$summary()

# primary parametersの名前を見る
names(fit_mcmc$EAP(pars = "parameters", drop = FALSE))

# MCMC drawの展開名を見る
dimnames(fit_mcmc$draws())[[3]]

# wrapperが作ったモデルコードを見る
print_code(mdl)

fit$estimate(pars = "parameters", drop = FALSE) は、推定されたprimary parametersだけをリストで返します。これを fixed や次のMCMCの init に渡すと、transformやgenerateで作られた派生量を誤って使うことを避けられます。

9.2 fixedの内部処理

BayesRTMBでは、ユーザーは自然尺度でパラメータを指定します。内部では、制約付きパラメータを非制約尺度へ変換し、必要なヤコビアン補正を推定経路に応じて処理します。

MCMCでは、非制約空間で事後分布をサンプリングするため、基本的にすべての制約変換に対応するヤコビアン補正が使われます。Laplace近似を含むMAPでは、random側または積分される側に必要な補正が使われます。

fixed でパラメータを固定した場合も、固定値のprior contributionや制約変換の扱いが内部で調整されます。そのため、fixed を使って作ったnested modelは、Bayes factorやMAP推定で一貫したrestricted modelとして扱えます。

9.3 fixedの制限

fixed はrestricted modelを作るための機能ですが、すべての部分固定が安全にできるわけではありません。

多変量priorや構造制約のあるparameterでは、一部の要素だけを固定すると、元の制約空間・ヤコビアン・prior contributionを正しく分解できないことがあります。BayesRTMBは危険なケースを検出した場合はエラーにします。Bayes factorを目的にする場合は、固定対象が明確なスカラーか、比較仮説に対応する縮小モデルとして自然に書ける形にするのが安全です。

9.4 optimize結果を固定値として使う

fit_map <- mdl$optimize(num_estimate = 4)

est <- fit_map$estimate(
  pars = "parameters",
  drop = FALSE
)

mdl_fixed <- mdl$fixed_model(
  fixed = list(
    sigma = est$sigma
  )
)

fit_fixed <- mdl_fixed$optimize()

特定の係数だけを固定する場合は、展開名を使います。

mdl_fixed <- mdl$fixed_model(
  fixed = list(
    "b[talk]" = 0,
    "b[perf]" = 0
  )
)

9.5 MCMC結果を固定値として使う

MCMCのEAPやMAPを固定値として使うこともできます。

fit_mcmc <- mdl$sample()

eap <- fit_mcmc$EAP(
  pars = "parameters",
  drop = FALSE
)

mdl_fixed <- fit_mcmc$model$fixed_model(
  fixed = list(
    sigma = eap$sigma
  )
)

fit_fixed <- mdl_fixed$sample()

ベクトルパラメータを丸ごと固定する場合は、そのパラメータ名に値を渡します。

eap <- fit_mcmc$EAP(pars = "parameters", drop = FALSE)

mdl_b_fixed <- fit_mcmc$model$fixed_model(
  fixed = list(
    b = eap$b
  )
)

一部の要素だけを固定する場合は、展開名を使います。

mdl_one_fixed <- fit_mcmc$model$fixed_model(
  fixed = list(
    "b[talk]" = 0
  )
)

9.6 fixedとBayes factor

bayes_factor() の fixed は、restricted modelを作るための便利なショートカットです。

bf <- fit_mcmc$bayes_factor(
  fixed = list("b[talk]" = 0),
  chains = 4,
  sampling = 4000
)

この指定は、次の処理をまとめて行うのと同じ考え方です。

mdl_null <- fit_mcmc$model$fixed_model(
  fixed = list("b[talk]" = 0)
)

fit_null <- mdl_null$sample(
  chains = 4,
  sampling = 4000
)

bf <- fit_mcmc$bayes_factor(
  comparison_fit = fit_null
)

10.1 setupに置くべき処理

データだけから決まる処理は、できるだけ setup に置きます。典型的には、次のような処理です。

• data frameの列を、モデル内で使う変数名に結び付ける。
• 欠測位置やindexを作る。
• デザイン行列、カテゴリ数、観測数などの定数を作る。
• AD対象外でよい補助関数を定義する。

setup = {
  Y <- response
  X <- model.matrix(~ group + score)

  obs <- which(!is.na(Y), arr.ind = TRUE)
  person_idx <- as.integer(obs[, "row"])
  item_idx <- as.integer(obs[, "col"])
  Y_obs <- Y[obs]
  N_obs <- length(Y_obs)

  center <- function(x) x - mean(x)
}

setup はADテープに記録されません。欠測位置、index、定数、デザイン行列の作成などに向いています。

10.2 transformに置くべき処理

パラメータに依存し、かつ model ブロックでも使う派生量は transform に置きます。典型例は、線形予測子、平均構造、確率、相関行列、標準化係数などです。

transform = {
  eta <- X %*% b + Intercept
  mu <- inv_logit(eta)
}

transform で作った量は、model ブロックの中でそのまま使えます。また、推定後にposterior drawごとに再計算できます。

fit$transformed_draws()
fit$estimate(pars = "transform")

transform 内で代入した変数は、デフォルトでは保存対象になります。途中計算を保存したくない場合は、report() で出力したい変数だけを明示します。

transform = {
  eta <- X %*% b + Intercept
  mu <- inv_logit(eta)
  report(mu)
}

この例では、eta は計算途中で使うだけなので保存されず、mu だけが transform 成分に保存されます。複数の量を保存したい場合は、report(mu, eta) のように指定します。

report() は「計算するかどうか」ではなく「fit objectに保存して、estimate()、draws()、transformed_draws() などから取り出せるようにするか」を制御するものです。report() しなかった途中変数も、そのブロック内や後続の model 計算には使えます。

10.3 generateに置くべき処理

推定後に必要だが、log posteriorの計算そのものには不要な量は generate に置きます。典型例は、posterior prediction、WAIC用のpointwise log_lik、予測確率、分類結果、回転後の量、解釈用の派生量などです。

generate = {
  log_lik <- normal_lpdf(Y, eta, sigma, sum = FALSE)
  y_rep <- rnorm(length(Y), eta, sigma)
  report(log_lik, y_rep)
}

generate で作った量は、model ブロックの尤度や事前分布には影響しません。推定後に generated_quantities() でdrawごとに計算され、estimate(pars = "generate") や draws(pars = ...) から取り出せます。

gq_code <- rtmb_code(
  generate = {
    log_lik <- normal_lpdf(Y, eta, sigma, sum = FALSE)
  }
)

fit$generated_quantities(gq_code)
fit$estimate(pars = "generate")
fit$draws(pars = "log_lik")

generate でも report() を使えます。report() を書かない場合は、ブロック内で代入した量が保存対象になります。保存したい量を絞る場合は、report(log_lik) のように明示します。

log_lik を保存しておくと WAIC() が使えます。WAIC用の log_lik は、全観測を足し合わせた1つの値ではなく、観測ごとのベクトルである必要があります。そのため normal_lpdf(..., sum = FALSE) のように書きます。

10.4 transformとgenerateの使い分け

迷った場合は、「この量がlikelihoodやpriorを計算するために必要か」で判断します。必要なら transform、推定後に見たいだけなら generate が基本です。

11.1 次元

3次元以上の配列も使えます。Dim(c(I, J, K)) のように、次元を整数ベクトルで指定します。

11.2 制約

制約付きパラメータは、内部では非制約尺度に変換されます。MCMCやLaplace近似では、必要なヤコビアン補正が推定経路に応じて使われます。

11.3 random = TRUE

random = TRUE は、そのパラメータをrandom effectとして扱う指定です。

parameters = {
  theta <- Dim(N_persons, random = TRUE)
  b <- Dim(N_items)
}

optimize(laplace = TRUE) では、random = TRUE のパラメータはLaplace近似で積分されます。sample() では、random effectも含めた事後サンプルを得られますが、metricの扱いではrandom側が特別に扱われることがあります。

12.1 連続分布

多くの単変量分布はベクトル化されています。

Y ~ normal(mu, sigma)

分布の引数は、名前の通りの尺度で指定します。たとえば lognormal(meanlog, sdlog) は、元の値ではなく、対数を取った値の平均と標準偏差を指定します。bernoulli_logit(eta) や categorical_logit(eta) は、確率ではなくlogit尺度の値を渡します。

観測ごとのlog likelihoodが必要な場合は、sum = FALSE を使います。

log_lik <- normal_lpdf(Y, mu, sigma, sum = FALSE)

sum = FALSE は、各観測のlog densityをベクトルとして返します。WAIC() やpointwise log likelihoodが必要なモデル比較では、この形で generate ブロックに保存します。

12.2 離散分布

離散分布のカテゴリ番号は、基本的にStanと同じく1始まりです。categorical()、ordered_logistic() などで0始まりのカテゴリを渡すと、意図しないモデルになります。

12.3 多変量・行列分布

12.4 混合分布・特殊分布

14.1 setupを活用する

パラメータに依存しない処理は setup に置きます。

setup = {
  obs <- which(!is.na(Y), arr.ind = TRUE)
  person_idx <- as.integer(obs[, "row"])
  item_idx <- as.integer(obs[, "col"])
  Y_obs <- Y[obs]
}

このような処理を model 内に置くと、ADテープ化のたびに不要な処理が記録され、遅くなります。

14.2 ベクトル化を優先する

同じ計算を書けるなら、明示的なループよりもベクトル化した行列演算を優先します。

transform = {
  eta <- X %*% b + Intercept
}

model = {
  Y ~ normal(eta, sigma)
}

ただし、ループ自体が禁止ではありません。AD対象の演算が微分可能であれば、for ループは使えます。

model = {
  for (i in 1:N_obs) {
    eta_i <- a[item_idx[i]] * (theta[person_idx[i]] - b[item_idx[i]])
    Y_obs[i] ~ bernoulli_logit(eta_i)
  }
}

14.3 apply系関数を避ける

apply()、sapply()、lapply() は、AD属性を落とすことがあります。rtmb_code() 内では、明示的なループかベクトル化を使うほうが安全です。

14.4 ifとifelseの注意

パラメータ値に依存する条件分岐は、微分可能性を壊すことがあります。

避けたい例:

model = {
  if (theta > 0) {
    y ~ normal(theta, 1)
  } else {
    y ~ normal(0, 1)
  }
}

必要な場合は、滑らかな近似、log_mix()、またはモデル構造の再定式化を検討します。

14.5 rtmb_vector()とrtmb_array()

代入でベクトルや配列を作る場合、numeric() や array() ではAD型がうまく保持されないことがあります。その場合は rtmb_vector() や rtmb_array() を使います。

transform = {
  eta <- rtmb_vector(0, length = N_obs)
  for (i in 1:N_obs) {
    eta[i] <- X[i, ] %*% b
  }
}

transform = {
  logit_x <- rtmb_array(0, dim = c(T, C, D))
  for (t in 1:T) {
    for (c in 1:C) {
      for (d in 1:D) {
        logit_x[t, c, d] <- alpha[d] + beta[c, d] * time[t]
      }
    }
  }
}

ただし、rtmb_vector() や rtmb_array() は代入を多用するモデルの補助です。行列演算やベクトル化で自然に書ける場合は、そちらを優先します。

14.6 テープ化時間とAD評価速度

テープ化時間は、モデルをRTMBのADオブジェクトとして記録する時間です。AD評価速度は、記録されたテープを使ってlog posteriorやgradientを繰り返し評価する速度です。

長いループや大量の代入は、テープ化を遅くすることがあります。行列演算にできる部分は行列演算にすると、テープ化時間とAD評価速度の両方が改善しやすくなります。

14.7 ADテープ内で避けたい書き方

モデルコードはRとして実行できても、ADテープとしては遅い、または不安定なことがあります。迷ったときは、次の表を確認します。

generate ブロックは推定後にdrawごとに評価する場所なので、posterior predictionの乱数生成や log_lik の保存に向いています。一方、model ブロックはlog posteriorそのものを定義する場所なので、同じパラメータなら同じ値を返す決定的な計算にします。

15.1 MCMCの結果から推定値だけ欲しい

fit <- mdl$sample()

fit$EAP(pars = "parameters")
fit$MAP(pars = "parameters")
fit$estimate(pars = "parameters", type = "EAP")

15.2 初期値をMAPから作ってMCMCを回したい

fit_map <- mdl$optimize(num_estimate = 4)

fit_mcmc <- mdl$sample(
  init = fit_map$estimate(pars = "parameters")
)

15.3 生成量だけを取り出したい

fit$generated_quantities({
  pred <- rnorm(length(Y), mu, sigma)
})

fit$estimate(pars = "generate")
fit$draws(pars = "pred")

15.4 random effectsを含めてdrawを取り出したい

draws_re <- fit$draws(
  inc_random = TRUE,
  inc_transform = FALSE,
  inc_generate = FALSE
)

15.5 VBの収束を見たい

fit_vb <- mdl$variational(iter = 10000, num_estimate = 4)

fit_vb$plot_elbo(tail_n = 1000)
fit_vb$diagnose()

ELBOが後半で水平になっていれば、VBの最適化は安定している可能性が高いです。ただし、VBは近似なので、重要な結論はMCMCで確認します。

15.6 classicで頻度主義的にモデル比較したい

fit0 <- rtmb_lm(sat ~ 1, data = debate, prior = prior_flat())$classic()
fit1 <- rtmb_lm(sat ~ talk + perf, data = debate, prior = prior_flat())$classic()

anova(fit0, fit1)
AIC(fit0, fit1)
BIC(fit0, fit1)

15.7 Bayes factorで係数の有無を比較したい

fit <- rtmb_lm(
  sat ~ talk * perf,
  data = debate,
  prior = prior_normal()
)$sample()

bf <- fit$bayes_factor(
  fixed = list("b[talk:perf]" = 0),
  chains = 4,
  sampling = 4000
)

15.8 fixed modelを明示的に作りたい

mdl_full <- rtmb_lm(
  sat ~ talk * perf,
  data = debate,
  prior = prior_normal()
)

mdl_null <- mdl_full$fixed_model(
  fixed = list("b[talk:perf]" = 0)
)

fit_full <- mdl_full$sample()
fit_null <- mdl_null$sample()

15.9 古いfit objectを新しいクラスに更新したい

fit_new <- upgrade_fit(fit_old)

保存済みfit objectを現在のクラス定義で使い直すための関数です。パッケージがバージョンアップした場合などに使います。

15.10 並列でMCMCを回したい

fit <- mdl$sample(
  parallel = TRUE,
  chains = 4,
  progress = "message"
)

parallel = TRUE にすると、chainを並列に実行します。progress = "message" は、Jobや一部の環境でも進捗を確認しやすい表示です。進捗表示が不要な場合は progress = "none" を使います。

通常は、モデルに必要なデータや関数を data と setup に入れておけば十分です。外部のオブジェクトや関数をどうしてもworkerへ渡したい場合は、推定時に globals = TRUE を指定します。

fit <- mdl$sample(
  parallel = TRUE,
  globals = TRUE
)