R Academy · Урок

Выборка MCMC и диагностика

Запускайте выборку, изучайте цепи и интерпретируйте диагностические показатели Rhat и ESS

«Выборка MCMC и диагностика» — бесплатный урок R Academy на CoddyKit. Это урок 3 из 4. Ты можешь прочитать весь урок бесплатно ниже — а потом практиковать его прямо в браузере с встроенным редактором кода и ИИ-репетитором 24/7. Это часть пути обучения R Academy, и твой прогресс синхронизируется между веб-версией и приложением CoddyKit. Курс R Academy содержит 4 уроков всего.

Что такое MCMC

Метод Монте-Карло по цепям Маркова (MCMC) — это семейство алгоритмов для получения выборок из распределения вероятностей, когда прямая выборка невозможна. В байесовской статистике MCMC используется для получения выборки из апостериорного распределения P(parameters | data).

RStan реализует метод выборки без разворота (NUTS) — современный алгоритм MCMC.

Простая модель 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 выборок после разогрева — всего 4000 выборок для 4 цепей.

# 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. Эти два диагностических показателя следует проверять в первую очередь.

# 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    1

Критерий сходимости Rhat

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) открывает интерактивное приложение Shiny с графиками трасс, апостериорными распределениями, графиками pairs и диагностикой NUTS в одном месте. Это наиболее полный инструмент для изучения подгонки 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
Увеличьте adapt_delta в направлении 1.0 (например, до 0.95) в control
Измените параметризацию — используйте нецентрированную параметризацию для иерархических моделей
Сделайте априорные распределения более узкими, если они слишком расплывчаты
Проверьте данные — выбросы или различия в масштабах вызывают проблемы у алгоритма выборки

# 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

Каков текущий стандартный порог Rhat, указывающий на сходимость модели Stan?

Повторение: выборка MCMC и диагностика

Основной рабочий процесс MCMC в RStan:

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/7) и разблокировать остальной курс R Academy, подпишись на CoddyKit PRO. Курс R Academy содержит 4 уроков всего.

Чему я научусь в уроке «Выборка MCMC и диагностика»?

Запускайте выборку, изучайте цепи и интерпретируйте диагностические показатели Rhat и ESS Ты практикуешь R Academy с помощью реального кода, который запускаешь прямо в браузере, и ИИ-репетитор 24/7 отвечает на твои вопросы во время урока.

Нужен ли мне опыт, чтобы начать R Academy?

Предыдущий опыт не требуется. R Academy на CoddyKit структурирован для всех уровней — от новичков до продвинутых, поэтому ты можешь начать отсюда или с самого начала и учиться в своем темпе. Это урок 3 из 4.

Сколько времени занимает урок «Выборка MCMC и диагностика»?

Большинство уроков CoddyKit занимают около 5–10 минут. Каждый из них компактный и интерактивный, поэтому ты постоянно делаешь прогресс и продолжаешь с того же места в веб-версии и приложении.

Можно ли писать и запускать код в этом уроке R Academy?

Да. Каждый урок R Academy включает встроенный редактор кода, поэтому ты пишешь и запускаешь реальный код прямо в браузере и получаешь моментальную обратную связь от AI — локальная установка не требуется.

Все уроки этого курса

← Назад к R Academy