RStan:MCMCにおけるパラメーター非中心化の効果

model_code <- "
data {
  int<lower=1> N;
  int<lower=1> J;
  vector[N] y;
  int<lower=1, upper=J> group_id[N];
}

parameters {
  real mu_theta;              // グループ平均の全体平均
  real<lower=0> sigma_theta;  // グループ平均の標準偏差
  vector[J] theta;            // 各グループの平均
  real<lower=0> sigma_y;      // 観測誤差の標準偏差
}

model {
  // 事前分布
  mu_theta ~ normal(0, 10);
  sigma_theta ~ cauchy(0, 2.5); // 弱情報事前分布
  sigma_y ~ cauchy(0, 2.5);

  // 階層構造
  theta ~ normal(mu_theta, sigma_theta);

  // 尤度
  y ~ normal(theta[group_id], sigma_y);
}
"

model_code_ncp <- "
data {
  int<lower=1> N;
  int<lower=1> J;
  vector[N] y;
  int<lower=1, upper=J> group_id[N];
}

parameters {
  real mu_theta;
  real<lower=0> sigma_theta;
  vector[J] theta_raw;        // 非中心化された「生の」パラメータ
  real<lower=0> sigma_y;
}

transformed parameters {
  vector[J] theta;
  // thetaをtheta_rawから変換して生成
  theta = mu_theta + theta_raw * sigma_theta;
}

model {
  // 事前分布
  mu_theta ~ normal(0, 10);
  sigma_theta ~ cauchy(0, 2.5);
  sigma_y ~ cauchy(0, 2.5);

  // 階層構造 (非中心化)
  theta_raw ~ normal(0, 1); // theta_rawが標準正規分布に従う

  // 尤度
  y ~ normal(theta[group_id], sigma_y);
}
"

# stanfit オブジェクトの読み込み
setwd(stan_output)
fit <- readRDS("stan_fit.rds")
fit_ncp <- readRDS("stan_fit_ncp.rds")

# 3. 診断情報の比較
# ------------------------------------
# 発散的移行の数を取得
sampler_params <- get_sampler_params(fit, inc_warmup = FALSE)
divergences <- sum(sapply(sampler_params, function(x) sum(x[, "divergent__"])))
cat("非中心化しないモデルの発散的移行の数:", divergences, "\n")

sampler_params_ncp <- get_sampler_params(fit_ncp, inc_warmup = FALSE)
divergences_ncp <- sum(sapply(sampler_params_ncp, function(x) sum(x[, "divergent__"])))
cat("非中心化モデルの発散的移行の数:", divergences_ncp, "\n")

# パラメータサマリーの比較（特にn_effとRhat）
summary <- summary(fit, pars = c("mu_theta", "sigma_theta", "sigma_y"))$summary
summary_ncp <- summary(fit_ncp, pars = c("mu_theta", "sigma_theta", "sigma_y"))$summary

cat("\n--- 非中心化しないモデルのサマリー ---\n")
print(summary)
cat("\n--- 非中心化モデルのサマリー ---\n")
print(summary_ncp)