RStan:MCMC収束の確認方法

1. 必要なパッケージの読み込み

rstanと、結果の可視化のためにbayesplotパッケージを使用します。

library(rstan)
library(bayesplot)

sapply(X = c("rstan", "bayesplot"), packageVersion)

stan_output <- "D:/stan_output"

$rstan
[1]  2 32  7

$bayesplot
[1]  1 13  0

2. サンプルデータの準備

平均10、標準偏差2の正規分布から100個のデータを生成します。このmu=10とsigma=2を、Rstanで推定するターゲットとします。

seed <- 20250626
set.seed(seed)
N <- 100
mu_true <- 10
sigma_true <- 2
y <- rnorm(N, mean = mu_true, sd = sigma_true)

1. Stanモデルの記述

観測データyが正規分布に従うと仮定し、その平均muと標準偏差sigmaを推定するモデルを記述します。

# Stanモデル
stan_model <- "
data {
  int<lower=0> N;  // データ数
  vector[N] y;     // 観測データ
}

parameters {
  real mu;               // 推定する平均
  real<lower=0> sigma;   // 推定する標準偏差（正の値のみ）
}

model {
  // 事前分布（弱い無情報事前分布）
  mu ~ normal(0, 100);
  sigma ~ cauchy(0, 5);

  // 尤度
  y ~ normal(mu, sigma);
}
"

2. MCMCサンプリングの実行

作成したモデルとデータを使って、MCMCサンプリングを実行します。ここでは4つのチェーンを走らせます。

rstan_options(auto_write = TRUE)
options(mc.cores = parallel::detectCores())

# Stanに渡すためのデータリストを作成
dat <- list(N = N, y = y)

# MCMCの実行
fit <- stan(
  model_code = stan_model, # Stanモデル
  data = dat, # データ
  iter = 2000, # サンプリング回数（ウォームアップ含む）
  warmup = 1000, # ウォームアップ（捨てる）期間
  chains = 4, # チェーンの数
  seed = seed # MCMCの再現性のための乱数シード
)

# stanfit オブジェクトの保存
setwd(stan_output)
saveRDS(object = fit, file = "stan_fit.rds")

• iter: 1チェーンあたりのサンプリング総数
• warmup: パラメータの初期値の影響をなくすための準備期間。この期間のサンプルは捨てられます。
• chains: 独立したMCMCをいくつ実行するか。複数（通常は3〜4以上）実行し、それらが同じ結果に収束するかを確認します。

1. R-hat (\(\\hat{R}\)​)

R-hatは、チェーン内分散とチェーン間分散を比較することで、チェーンが同じ分布に収束したかを評価する指標です。

• 判断基準: Rhat < 1.1（より厳しく 1.05 や 1.01 を基準にすることもあります）
• 理想: Rhatが1.0に近いほど良い。

今回の結果の解釈: 上の結果を見ると、muとsigmaのRhatは両方とも1です。これは、4つのチェーンが同じ分布にうまく収束したことを示しており、収束は良好と判断できます。

2. 有効サンプルサイズ (n_eff)

n_eff (Effective Sample Size)は、MCMCサンプル間の自己相関を考慮した、実質的に独立なサンプルの数を示します。

• 判断基準: n_effが小さすぎると（例えば100未満）、事後分布の推定が不安定になる可能性があります。多ければ多いほど良いです。
• 自己相関との関係: MCMCサンプルは通常、自己相関を持ちます（1つ前のサンプルと似ている）。自己相関が高いとn_effは小さくなります。

今回の結果の解釈: n_effはmuで3150 、sigmaで3407です。総サンプル数4000に対して十分に大きな値であり、事後分布を要約するのに十分なサンプルが得られていると言えます。

事後分布のプロット

各チェーンから得られた事後分布を重ねて描画し、形状が一致しているかを確認します。

# ヒストグラムで事後分布を確認
mcmc_hist(fit, pars = c("mu", "sigma"))

# 密度プロットで事後分布を確認
mcmc_dens_overlay(fit, pars = c("mu", "sigma"))

解釈: 4つのチェーンから得られた分布がほぼ重なっており、Rhatが1であったことを視覚的にも裏付けています。

自己相関プロット

サンプルの自己相関がどのくらいの速さで減衰するかを確認します。減衰が速い（すぐに0に近づく）ほど、効率的なサンプリングができており、n_effが大きくなります。

mcmc_acf(fit, pars = c("mu", "sigma")) + ggplot2::scale_x_continuous(breaks = scales::pretty_breaks(10))

解釈: ラグ1以降の自己相関が急速に0に近づいています。これはサンプルの自己相関が低く、効率的なサンプリングが行われたことを示します。