R Academy · Урок

Проверка апостериорных предсказаний

Проверяйте соответствие модели, сравнивая распределения смоделированных и наблюдаемых данных

Урок 4 из 413 шагов

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

Что такое апостериорные прогнозные проверки

После подгонки байесовской модели необходимо спросить: генерирует ли эта модель данные, похожие на наблюдавшиеся данные? Апостериорные прогнозные проверки (PPC) отвечают на этот вопрос, моделируя повторные наборы данных yrep из апостериорного распределения и сравнивая их графически с наблюдаемыми данными y.

Извлечение апостериорных выборок

extract(fit, pars='mu') возвращает именованный список; $mu — числовой вектор всех выборок этого параметра после разогрева. При 4 цепях и 1000 итераций после разогрева в каждой цепи получается 4000 выборок.

# library(rstan)
# mu_samples    <- extract(fit, pars = 'mu')$mu
# sigma_samples <- extract(fit, pars = 'sigma')$sigma
#
# cat('Samples drawn:', length(mu_samples), '
')
# cat('Posterior mean of mu:', mean(mu_samples), '
')
# cat('90% CI:', quantile(mu_samples, c(0.05, 0.95)), '
')

Генерация yrep из апостериорных выборок

Для каждого апостериорного извлечения (mu_s, sigma_s) смоделируйте повторный набор данных того же размера, что и исходные данные. Сохраните эти наборы в матрице yrep, где каждая строка представляет один смоделированный набор данных.

# y <- c(2.1, 1.8, 2.4, 1.9, 2.3, 2.0, 1.7, 2.2)
# n_obs <- length(y)
# S     <- length(mu_samples)  # 4000 posterior draws
#
# yrep <- matrix(NA, nrow = S, ncol = n_obs)
# for (s in seq_len(S)) {
#   yrep[s, ] <- rnorm(n_obs, mean = mu_samples[s], sd = sigma_samples[s])
# }
# dim(yrep)  # [4000, 8]

ppc_dens_overlay() — сравнение плотностей

bayesplot::ppc_dens_overlay(y, yrep[1:50,]) накладывает ядерную плотность наблюдаемых данных (тёмная линия) на плотности 50 случайно выбранных смоделированных наборов данных (светлые линии). Хорошее соответствие модели означает, что тёмная линия находится внутри облака светлых линий.

# library(bayesplot)
#
# ppc_dens_overlay(y, yrep[1:50, ])
#
# Interpretation:
# - Dark line (y_obs) surrounded by light lines (yrep): good fit
# - Dark line systematically outside the cloud: model misfit
# - Light lines much wider than dark: overdispersed model
# - Light lines much narrower than dark: underdispersed model

ppc_stat() — проверка статистики

ppc_stat(y, yrep, stat = 'mean') показывает гистограмму проверочной статистики (например, среднего), вычисленной для каждого смоделированного набора данных, с вертикальной линией на уровне наблюдаемой статистики. Если наблюдаемое значение попадает в основную часть гистограммы, модель хорошо описывает этот аспект данных.

# library(bayesplot)
#
# ppc_stat(y, yrep, stat = 'mean')   # does model capture the mean?
# ppc_stat(y, yrep, stat = 'sd')     # does model capture spread?
# ppc_stat(y, yrep, stat = 'max')    # does model capture extremes?
#
# If observed stat is in the tail of the histogram,
# the model fails to reproduce that statistic.

Байесовское p-значение

Байесовское p-значение (апостериорное прогнозное p-значение) — это доля смоделированных наборов данных, в которых проверочная статистика оказывается более экстремальной, чем наблюдаемое значение. Значения около 0.5 указывают на хорошую калибровку; значения около 0 или 1 говорят о несоответствии модели по этой статистике.

# Bayesian p-value for the mean:
# obs_mean <- mean(y)
# rep_means <- apply(yrep, 1, mean)
# pval <- mean(rep_means >= obs_mean)
# cat('Bayesian p-value (mean):', round(pval, 3), '
')
# # 0.5 is perfect; < 0.05 or > 0.95 suggests misfit

Другие функции PPC в bayesplot

bayesplot предлагает множество визуализаций PPC помимо наложения плотностей:

  • ppc_hist(y, yrep[1:8,]) — сетка гистограмм
  • ppc_scatter_avg(y, yrep) — диаграмма рассеяния наблюдаемых значений и среднего yrep
  • ppc_intervals(y, yrep) — интервалы неопределённости вокруг каждого наблюдения
  • ppc_rootogram(y, yrep) — для данных-счётчиков
# library(bayesplot)
#
# # Grid of 8 simulated histograms vs the observed
# ppc_hist(y, yrep[1:8, ])
#
# # Scatter: y_obs (x) vs mean of yrep (y) — should hug diagonal
# ppc_scatter_avg(y, yrep)
#
# # 50% and 90% posterior predictive intervals around each y_i
# ppc_intervals(y, yrep)

Интерпретация графиков PPC — несоответствие модели

Распространённые признаки несоответствия и их причины:

  • yrep слишком широко распределён — априорное распределение слишком расплывчато или модель обладает избыточной дисперсией
  • yrep смещён — выбрано неправильное семейство распределений правдоподобия (например, нормальное для асимметричных данных)
  • yrep не воспроизводит мультимодальность — необходима смесь распределений
  • yrep не воспроизводит экстремальные значения — необходимо распределение с тяжёлыми хвостами
# Example: if data has a long right tail but yrep does not,
# consider switching:
# y ~ normal(mu, sigma)  =>  y ~ student_t(nu, mu, sigma)
#
# Or for count data:
# y ~ poisson(lambda)  =>  y ~ neg_binomial_2(mu, phi)  (overdispersion)
cat('PPCs guide model improvement by revealing specific failure modes
')

PPC в блоке модели Stan

Можно генерировать yrep непосредственно в Stan с помощью блока generated quantities. Это избавляет от повторного извлечения параметров в R и эквивалентно ему с вычислительной точки зрения.

# Stan model with generated quantities:
# '
# generated quantities {
#   array[N] real y_rep;
#   for (n in 1:N) {
#     y_rep[n] = normal_rng(mu, sigma);
#   }
# }
# '
# Then extract in R:
# yrep <- extract(fit, pars = 'y_rep')$y_rep  # [S, N] matrix

Перекрёстная проверка с исключением одного наблюдения

Помимо PPC, loo::loo(fit) вычисляет перекрёстную проверку с исключением одного наблюдения для сравнения конкурирующих моделей. Предпочтение отдаётся модели с более высоким ELPD (ожидаемой логарифмической плотностью прогнозирования). Используйте loo::loo_compare(loo1, loo2), чтобы ранжировать модели.

# library(loo)
# loo1 <- loo(fit1)  # normal model
# loo2 <- loo(fit2)  # student-t model
#
# comparison <- loo_compare(loo1, loo2)
# print(comparison)
#
# Model with elpd_diff > 0 is preferred
# se_diff > |elpd_diff| means difference is not reliable

Рекомендации по применению PPC

Следуйте этим рекомендациям для строгой апостериорной прогнозной проверки:

  • Всегда начинайте с ppc_dens_overlay() как с общей проверки здравого смысла
  • Затем используйте предметно-ориентированные статистики (ppc_stat()), соответствующие целям вашего анализа
  • Для визуальных проверок используйте не менее 50 извлечений yrep, а для байесовских p-значений — все 4000
  • Неудачные PPC помогают улучшать модель — это диагностический результат, а не ошибка
# Workflow:
# 1. Fit model -> extract() -> generate yrep matrix
# 2. ppc_dens_overlay(y, yrep[1:50,])  -- visual global check
# 3. ppc_stat(y, yrep, stat='mean')    -- check mean
# 4. ppc_stat(y, yrep, stat='sd')      -- check spread
# 5. ppc_stat(y, yrep, stat='max')     -- check tails
# 6. If misfit found -> revise model -> refit -> re-check

Быстрая проверка: интерпретация байесовского p-значения

Что означает байесовское p-значение 0.03 для статистики максимума с точки зрения модели?

Повторение: апостериорные прогнозные проверки

Рабочий процесс PPC в RStan и bayesplot:

  • Извлеките выборки: extract(fit, pars='mu')$mu
  • Сгенерируйте yrep: переберите апостериорные извлечения, вызывая rnorm(n, mu_s, sigma_s)
  • Общая проверка: ppc_dens_overlay(y, yrep[1:50,])
  • Проверка статистики: ppc_stat(y, yrep, stat='mean')
  • Байесовское p-значение: mean(apply(yrep,1,stat) >= stat(y)) — значение около 0.5 является хорошим
  • Для сравнения моделей используйте loo::loo_compare()
Можно начать бесплатно

Изучай R с ИИ-репетитором — бесплатно

Пиши и запускай код прямо в браузере, получай мгновенную помощь от ИИ-репетитора 24/7 и продолжи учиться на сайте или в приложении.

Курсы
43
Уроки
159

Часто задаваемые вопросы

Урок «Проверка апостериорных предсказаний» бесплатный?

Да — полный текст урока «Проверка апостериорных предсказаний» бесплатно доступен здесь в веб-версии. Чтобы практиковать его интерактивно (встроенный редактор кода и ИИ-репетитор 24/7) и разблокировать остальной курс R Academy, подпишись на CoddyKit PRO. Курс R Academy содержит 4 уроков всего.

Чему я научусь в уроке «Проверка апостериорных предсказаний»?

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

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

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

Сколько времени занимает урок «Проверка апостериорных предсказаний»?

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

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

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

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

  1. Введение в байесовское мышление
  2. Написание моделей Stan в R
  3. Выборка MCMC и диагностика
  4. Проверка апостериорных предсказаний
← Назад к R Academy