Проверка апостериорных предсказаний
Проверяйте соответствие модели, сравнивая распределения смоделированных и наблюдаемых данных
«Проверка апостериорных предсказаний» — бесплатный урок 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 modelppc_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)— диаграмма рассеяния наблюдаемых значений и среднего yrepppc_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 — локальная установка не требуется.
Все уроки этого курса
- Введение в байесовское мышление
- Написание моделей Stan в R
- Выборка MCMC и диагностика
- Проверка апостериорных предсказаний