R Academy · Lekcja

Próbkowanie MCMC i diagnostyka

Uruchamiaj próbkowanie, analizuj łańcuchy i interpretuj diagnostyki Rhat oraz ESS.

Lekcja 3 z 413 kroki

Próbkowanie MCMC i diagnostyka to bezpłatna lekcja R Academy na CoddyKit. To lekcja 3 z 4. Możesz przeczytać całą lekcję poniżej za darmo — a potem ćwiczyć ją interaktywnie w przeglądarce z wbudowanym edytorem kodu i tutorem AI dostępnym 24/7. To część ścieżki edukacyjnej R Academy, a Twój postęp synchronizuje się między webem a aplikacją CoddyKit. Kurs R Academy zawiera 4 lekcji w sumie.

Czym jest MCMC?

Monte Carlo z łańcuchem Markowa (MCMC) to rodzina algorytmów służących do losowania próbek z rozkładu prawdopodobieństwa, gdy bezpośrednie próbkowanie jest niemożliwe. W statystyce bayesowskiej MCMC losuje próbki z rozkładu a posteriori P(parameters | data).

RStan implementuje No-U-Turn Sampler (NUTS), najnowocześniejszy algorytm MCMC.

Prosty model Stan

Model Stan to blok tekstu definiujący typy danych, parametry i logarytm rozkładu a posteriori. Najprostszy model estymuje średnią rozkładu normalnego o znanej wariancji.

# 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
# }
# '

Wywoływanie stan() w celu próbkowania

stan(model_code=..., data=..., chains=4, iter=2000, warmup=1000) kompiluje model (jeden raz), a następnie losuje próbki. Przy wartościach iter=2000 i warmup=1000 każdy łańcuch generuje 1000 próbek po fazie rozgrzewki — łącznie 4000 próbek z 4 łańcuchów.

# 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) — tabela podsumowania

print(fit) wyświetla tabelę podsumowania dla każdego parametru, zawierającą średnią a posteriori, odchylenie standardowe, kwantyle, Rhat i n_eff. Te dwa diagnostyczne wskaźniki należy sprawdzić w pierwszej kolejności.

# 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

Kryterium zbieżności Rhat

Rhat (współczynnik redukcji potencjalnej skali) porównuje wariancję wewnątrz łańcuchów z wariancją między łańcuchami. Wartości bliskie 1,0 wskazują, że wszystkie łańcuchy zbiegały do tego samego rozkładu.

Rhat < 1.01 — zbieżność osiągnięta (obecny standard)
Rhat > 1.01 — łańcuchy nie mieszają się; należy wykonać więcej iteracji
Rhat > 1.1 — poważny problem ze zbieżnością

# 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 — efektywna liczebność próby

n_eff (efektywna liczebność próby) uwzględnia autokorelację między kolejnymi próbkami MCMC. Próbki skorelowane zawierają mniej informacji niż próbki niezależne.

n_eff bliskie łącznej liczbie iteracji — próbki prawie niezależne, wynik doskonały
n_eff / total_samples > 0.1 — ogólnie akceptowalne
Bardzo niskie n_eff — silna autokorelacja; należy rozważyć zmianę parametryzacji modelu

# 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() do wizualnej oceny zbieżności

traceplot(fit, pars = 'mu') przedstawia wartości próbkowane dla mu w kolejnych iteracjach każdego łańcucha. Łańcuchy, które osiągnęły zbieżność, przypominają rozmytą gąsienicę — wszystkie się nakładają, bez trendów ani dryfu. Łańcuchy rozbieżne błądzą lub pozostają od siebie oddzielone.

# 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() do oceny korelacji a posteriori

pairs(fit, pars = c('mu', 'sigma')) wyświetla macierz wykresów rozrzutu próbek z rozkładu a posteriori. Ujawnia korelacje między parametrami i wyróżnia rozbieżne przejścia (zaznaczone na czerwono), które wskazują obszary sprawiające próbkowaczowi trudność.

# 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!
')

Wyodrębnianie próbek a posteriori

extract(fit, pars = 'mu')$mu zwraca wektor liczbowy wszystkich próbek parametru mu po fazie rozgrzewki. Użyj tych próbek, aby obliczyć dowolne podsumowanie rozkładu a posteriori: średnią, przedziały wiarygodności lub prawdopodobieństwo warunku.

# 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')

Uruchamianie ShinyStan do interaktywnej diagnostyki

shinystan::launch_shinystan(fit) otwiera interaktywną aplikację Shiny zawierającą wykresy śladu, rozkłady a posteriori, wykresy pairs i diagnostykę NUTS w jednym miejscu. To najbardziej wszechstronne narzędzie do badania dopasowania 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

Typowe sposoby poprawy zbieżności

Gdy Rhat > 1.01 lub widoczne są rozbieżności:

Zwiększ wartości iter i warmup
Zwiększ adapt_delta w kierunku 1,0 (np. do 0,95) w obiekcie control
Zmień parametryzację — użyj parametryzacji niecentrowanej dla modeli hierarchicznych
Zaostrz rozkłady a priori, jeśli są zbyt rozproszone
Sprawdź dane — wartości odstające lub różnice skal powodują problemy próbkowacza

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

Szybkie sprawdzenie: próg Rhat

Jaki jest obecny standardowy próg Rhat wskazujący, że model Stan osiągnął zbieżność?

Odpowiedz w interaktywnej lekcji

Podsumowanie: próbkowanie MCMC i diagnostyka

Przebieg pracy z MCMC w RStan:

stan(model_code=..., data=..., chains=4, iter=2000, warmup=1000) dopasowuje model
print(fit) wyświetla Rhat i n_eff — podstawowe wskaźniki diagnostyki zbieżności
Rhat < 1.01 i n_eff / total > 0.1 wskazują na poprawnie zachowującą się próbę
traceplot() — wizualna kontrola mieszania; pairs() — ujawnia problemy z geometrią rozkładu a posteriori
extract(fit, pars='mu')$mu — dostęp do surowych próbek a posteriori
shinystan::launch_shinystan(fit) — kompleksowa interaktywna diagnostyka

Bezpłatny start

Ucz się R dzięki korepetycjom AI — za darmo

Pisz i uruchamiaj kod w przeglądarce, otrzymuj natychmiastową pomoc od korepetytora AI dostępnego 24/7 i kontynuuj naukę w sieci lub w aplikacji.

Kursy: 43
Lekcje: 159

Zacznij naukę bezpłatnie Wyświetl pełny kurs

Często zadawane pytania

Czy lekcja „Próbkowanie MCMC i diagnostyka” jest bezpłatna?

Tak — pełny tekst „Próbkowanie MCMC i diagnostyka” jest dostępny za darmo tutaj w sieci. Aby ćwiczyć ją interaktywnie (wbudowany edytor kodu i tutor AI dostępny 24/7) i odblokować resztę kursu R Academy, przejdź na CoddyKit PRO. Kurs R Academy zawiera 4 lekcji w sumie.

Co nauczysz się w „Próbkowanie MCMC i diagnostyka”?

Uruchamiaj próbkowanie, analizuj łańcuchy i interpretuj diagnostyki Rhat oraz ESS. Ćwiczysz R Academy z praktycznym kodem, który uruchamiasz bezpośrednio w przeglądarce, a tutor AI dostępny 24/7 odpowiada na Twoje pytania podczas pracy nad lekcją.

Czy potrzebuję doświadczenia, aby zacząć R Academy?

Nie wymagamy żadnego doświadczenia. R Academy w CoddyKit jest strukturyzowany dla początkujących i zaawansowanych użytkowników, więc możesz zacząć tutaj lub od początku i uczyć się w swoim tempie. To lekcja 3 z 4.

Ile czasu zajmuje lekcja „Próbkowanie MCMC i diagnostyka”?

Większość lekcji CoddyKit trwa około 5–10 minut. Każda lekcja to mały, interaktywny krok, dzięki czemu robisz systematyczne postępy i zawsze wracasz dokładnie do tego samego miejsca — na webie i w aplikacji.

Czy mogę pisać i uruchamiać kod w tej lekcji R Academy?

Tak. Każda lekcja R Academy zawiera wbudowany edytor kodu, więc piszesz i uruchamiasz prawdziwy kod bezpośrednio w przeglądarce i od razu otrzymujesz sprzężenie zwrotne od AI — bez konfiguracji na komputerze.

Wszystkie lekcje w tym kursie

← Powrót do R Academy