MCMC サンプリングと診断
サンプリングを実行し、チェーンを確認して、Rhat と ESS の診断結果を解釈します。
「MCMC サンプリングと診断」はCoddyKit上の無料R Academyレッスンです。 これはレッスン3/4です。 下記で完全なレッスンを無料で読むことができます。その後、ブラウザ内の組み込みコードエディタと24時間対応のAIチューターでハンズオン演習できます。 これはR Academy学習パスの一部であり、ウェブとCoddyKitアプリ全体で進捗が同期されます。 R Academyコースには全4レッスンが含まれています。
MCMCとは
マルコフ連鎖モンテカルロ法（MCMC）は、直接サンプリングができない確率分布からサンプルを抽出するためのアルゴリズム群です。ベイズ統計では、事後分布 P(parameters | data)からサンプルを抽出します。
RStanは、最先端のMCMCアルゴリズムであるNo-U-Turn Sampler（NUTS）を実装しています。
単純なStanモデル
Stanモデルは、データ型、パラメータ、対数事後分布を定義するテキストブロックです。最も単純なモデルでは、分散が既知の正規分布の平均を推定します。
# library(rstan)
#
# stan_code <- '
# data {
# int<lower=0> N;
# vector[N] y;
# }
# parameters {
# real mu;
# real<lower=0> sigma;
# }
# model {
# mu ~ normal(0, 10); // prior
# sigma ~ exponential(1); // prior
# y ~ normal(mu, sigma); // likelihood
# }
# 'stan()を呼び出してサンプリングする
stan(model_code=..., data=..., chains=4, iter=2000, warmup=1000)はモデルを一度コンパイルしてから、サンプルを抽出します。iter=2000、warmup=1000の場合、各チェーンからウォームアップ後のサンプルが1000個生成され、4チェーン全体では4000個になります。
# library(rstan)
# options(mc.cores = parallel::detectCores())
#
# y <- c(2.1, 1.8, 2.4, 1.9, 2.3, 2.0, 1.7, 2.2)
# stan_data <- list(N = length(y), y = y)
#
# fit <- stan(
# model_code = stan_code,
# data = stan_data,
# chains = 4,
# iter = 2000,
# warmup = 1000,
# seed = 42
# )print(fit) — 要約テーブル
print(fit)は、各パラメータについて、事後平均、標準偏差、分位点、Rhat、n_effを含む要約テーブルを表示します。この2つの診断指標を最初に確認してください。
# print(fit)
#
# Example output:
# mean se_mean sd 2.5% 25% 50% 75% 97.5% n_eff Rhat
# mu 2.05 0.00 0.15 1.76 1.95 2.05 2.15 2.34 3842 1
# sigma 0.22 0.00 0.06 0.13 0.18 0.21 0.25 0.37 3521 1
# lp__ 4.38 0.02 1.01 1.60 3.91 4.71 5.18 5.50 2148 1Rhatによる収束判定
Rhat（潜在スケール縮小係数）は、チェーン内の分散とチェーン間の分散を比較します。値が1.0に近ければ、すべてのチェーンが同じ分布に収束したことを示します。
- Rhat < 1.01 — 収束済み（現在の基準）
- Rhat > 1.01 — チェーンが混ざっていません。反復回数を増やしてください
- Rhat > 1.1 — 深刻な収束問題があります
# Check Rhat for all parameters:
# s <- summary(fit)$summary
# rhat_vals <- s[, 'Rhat']
# cat('Max Rhat:', max(rhat_vals, na.rm = TRUE), '
')
# if (any(rhat_vals > 1.01, na.rm = TRUE)) {
# warning('Convergence issue detected!')
# } else {
# cat('All Rhat < 1.01 — chains converged
')
# }n_eff — 有効サンプルサイズ
n_eff（有効サンプルサイズ）は、MCMCサンプル間の自己相関を考慮した値です。相関のあるサンプルは、独立したサンプルよりも含む情報量が少なくなります。
- n_effが総反復回数に近い — ほぼ独立したサンプルで、非常に良好です
- n_eff / total_samples > 0.1 — 一般に許容範囲です
- n_effが非常に小さい — 自己相関が強いため、モデルの再パラメータ化を検討してください
# s <- summary(fit)$summary
# n_eff_vals <- s[, 'n_eff']
# total_samples <- 4 * 1000 # chains * post-warmup iter
# ratio <- n_eff_vals / total_samples
# cat('n_eff ratio (mu) :', round(ratio['mu'], 2), '
')
# cat('n_eff ratio (sigma):', round(ratio['sigma'], 2), '
')traceplot()による視覚的な収束確認
traceplot(fit, pars = 'mu')は、各チェーンについて、反復ごとのmuのサンプル値をプロットします。収束したチェーンは、ぼんやりした毛虫のように見えます。つまり、トレンドやドリフトがなく、すべてのチェーンが重なります。収束していないチェーンは、さまよったり、互いに離れたままになったりします。
# library(rstan)
#
# traceplot(fit, pars = c('mu', 'sigma'), inc_warmup = FALSE)
#
# Good traceplot characteristics:
# - All 4 chains overlapping completely (same range)
# - No visible drift or trend
# - Rapid mixing (values jump around quickly)
# - No flat regions (stuck sampler)
cat('A healthy traceplot looks like a fuzzy caterpillar
')pairs()による事後分布の相関確認
pairs(fit, pars = c('mu', 'sigma'))は、事後サンプルの散布図行列を表示します。パラメータ間の相関を明らかにし、サンプラーが苦戦している領域を示す発散遷移（赤で表示）を見つけるのに役立ちます。
# pairs(fit, pars = c('mu', 'sigma'))
#
# What to look for:
# - Elliptical clouds: mild correlation (OK)
# - Banana / funnel shapes: reparameterization needed
# - Red dots (divergences): geometry problem in posterior
# => increase adapt_delta: stan(..., control=list(adapt_delta=0.95))
cat('Red dots in pairs() indicate divergent transitions — investigate!
')事後サンプルの抽出
extract(fit, pars = 'mu')$muは、muのウォームアップ後の全サンプルを含む数値ベクトルを返します。これらのサンプルを使って、平均、信用区間、条件を満たす確率など、任意の事後分布の要約を計算できます。
# mu_samples <- extract(fit, pars = 'mu')$mu
# cat('Posterior mean :', mean(mu_samples), '
')
# cat('95% CI:', quantile(mu_samples, c(0.025, 0.975)), '
')
# cat('P(mu > 2):', mean(mu_samples > 2), '
')
# hist(mu_samples, main = 'Posterior of mu', xlab = 'mu', col = 'steelblue')ShinyStanを起動して対話的に診断する
shinystan::launch_shinystan(fit)は、traceplot、事後分布、pairsプロット、NUTS診断を1か所で確認できる対話的なShinyアプリを開きます。RStanのフィット結果を調べるための、最も包括的なツールです。
# install.packages('shinystan')
# library(shinystan)
#
# shinystan::launch_shinystan(fit)
#
# ShinyStan tabs:
# - Diagnose: Rhat, n_eff, divergences, energy
# - Explore: marginal posteriors, scatter plots
# - Model: Stan code, data
# - NUTS: step size, tree depth per chain収束問題への一般的な対処法
Rhat > 1.01の場合、または発散が見られる場合：
iterとwarmupを増やしますcontrolでadapt_deltaを1.0に近づけます（例：0.95）- 再パラメータ化します — 階層モデルでは非中心化パラメータ化を使用します
- 事前分布が広すぎる場合は、より強く設定します
- データを確認します — 外れ値やスケールの違いがサンプラーの問題を引き起こすことがあります
# Re-run with higher adapt_delta to reduce divergences:
# fit2 <- stan(
# model_code = stan_code,
# data = stan_data,
# chains = 4,
# iter = 4000,
# warmup = 2000,
# control = list(adapt_delta = 0.95, max_treedepth = 12),
# seed = 42
# )確認問題：Rhatのしきい値
Stanモデルが収束したことを示す、Rhatの現在の標準的なしきい値はいくつですか。
MCMCサンプリングと診断のまとめ
RStanにおけるMCMCの基本ワークフロー：
stan(model_code=..., data=..., chains=4, iter=2000, warmup=1000)でモデルをフィッティングしますprint(fit)でRhatとn_effを確認します — 主要な収束診断指標です- Rhat < 1.01かつn_eff / total > 0.1なら、サンプルは概ね良好です
traceplot()— 混合状態を視覚的に確認します。pairs()— 事後分布の幾何学的な問題を明らかにしますextract(fit, pars='mu')$mu— 生の事後サンプルにアクセスしますshinystan::launch_shinystan(fit)— 包括的な対話型診断を行います
よくある質問
「MCMC サンプリングと診断」レッスンは無料ですか？
はい。「MCMC サンプリングと診断」の完全なテキストはこのウェブで無料で読めます。インタラクティブに演習し（組み込みコードエディタと24時間対応のAIチューター）、R Academyコースの残りをアンロックするには、CoddyKit PROにアップグレードしてください。 R Academyコースには全4レッスンが含まれています。
「MCMC サンプリングと診断」で何を学びますか？
サンプリングを実行し、チェーンを確認して、Rhat と ESS の診断結果を解釈します。 ブラウザで直接実行するハンズオンコードでR Academyを演習し、24時間対応のAIチューターがレッスンを進める中での質問に答えます。
R Academyを始めるのに経験は必要ですか？
事前経験は必要ありません。CoddyKitのR Academyは初級者から上級者向けに構成されているため、ここから始めるか最初から始めて、自分のペースで進むことができます。 これはレッスン3/4です。
「MCMC サンプリングと診断」レッスンにはどのくらい時間がかかりますか？
ほとんどのCoddyKitレッスンは約5～10分かかります。各レッスンはコンパクトでインタラクティブなので、着実に進歩し、ウェブとアプリ全体で正確に前回の場所から再開できます。
このR Academyレッスンでコードを書いて実行できますか？
はい。すべてのR Academyレッスンに組み込みコードエディタが含まれているため、ブラウザでリアルコードを書いて実行し、即座のAIフィードバックを取得できます。ローカル設定は不要です。
このコースのすべてのレッスン
- ベイズ的思考の基礎
- R で Stan モデルを書く
- MCMC サンプリングと診断
- 事後予測チェック