R Academy · Lezione

Campionamento MCMC e diagnostica

Esegua il campionamento, esamini le catene e interpreti le diagnostiche Rhat ed ESS

Campionamento MCMC e diagnostica è una lezione R Academy gratuita su CoddyKit. Questa è la lezione 3 di 4. Puoi leggere la lezione completa qui gratuitamente — poi esercitati direttamente nel browser con un editor di codice integrato e un tutor IA disponibile 24/7. Fa parte del percorso di apprendimento R Academy, e i tuoi progressi si sincronizzano tra il web e l'app CoddyKit. Il corso R Academy include 4 lezioni in totale.

Che cos'è MCMC?

Markov Chain Monte Carlo (MCMC) è una famiglia di algoritmi per estrarre campioni da una distribuzione di probabilità quando il campionamento diretto è impossibile. Nella statistica bayesiana, MCMC estrae campioni dalla distribuzione a posteriori P(parameters | data).

RStan implementa il No-U-Turn Sampler (NUTS), un algoritmo MCMC all'avanguardia.

Un semplice modello Stan

Un modello Stan è un blocco di testo che definisce i tipi di dati, i parametri e il logaritmo della distribuzione a posteriori. Il modello più semplice stima la media di una distribuzione normale con varianza nota.

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

Chiamare stan() per il campionamento

stan(model_code=..., data=..., chains=4, iter=2000, warmup=1000) compila il modello (una sola volta), quindi estrae i campioni. Con iter=2000 e warmup=1000, ogni catena produce 1000 campioni dopo il warmup, per un totale di 4000 campioni nelle 4 catene.

# 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) — La tabella riepilogativa

print(fit) mostra una tabella riepilogativa per ciascun parametro, con media a posteriori, deviazione standard, quantili, Rhat e n_eff. Questi due diagnostici sono i primi elementi da controllare.

# 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

Il criterio di convergenza Rhat

Rhat (fattore di riduzione della scala potenziale) confronta la varianza all'interno delle catene con la varianza tra le catene. Valori vicini a 1.0 indicano che tutte le catene sono arrivate alla stessa distribuzione.

Rhat < 1.01 — convergenza raggiunta (standard attuale)
Rhat > 1.01 — le catene non hanno effettuato il mixing; esegua più iterazioni
Rhat > 1.1 — grave problema di convergenza

# 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 — Dimensione effettiva del campione

n_eff (dimensione effettiva del campione) tiene conto dell'autocorrelazione tra campioni MCMC successivi. I campioni correlati contengono meno informazioni rispetto a quelli indipendenti.

n_eff vicino al numero totale di iterazioni — campioni quasi indipendenti, eccellente
n_eff / total_samples > 0.1 — generalmente accettabile
n_eff molto basso — autocorrelazione elevata; consideri di riparametrizzare il modello

# 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() per la convergenza visiva

traceplot(fit, pars = 'mu') traccia i valori campionati di mu nelle varie iterazioni per ciascuna catena. Le catene convergenti assomigliano a un bruco sfocato: tutte le catene si sovrappongono, senza tendenze o derive. Le catene divergenti si spostano senza una direzione stabile o rimangono separate.

# 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() per le correlazioni a posteriori

pairs(fit, pars = c('mu', 'sigma')) mostra una matrice di grafici a dispersione dei campioni a posteriori. Rivela le correlazioni tra i parametri e mette in evidenza le transizioni divergenti (tracciate in rosso), che indicano le regioni con cui il sampler ha difficoltà.

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

Estrarre i campioni a posteriori

extract(fit, pars = 'mu')$mu restituisce un vettore numerico di tutti i campioni di mu successivi al warmup. Usi questi campioni per calcolare qualsiasi riepilogo a posteriori: media, intervalli di credibilità o probabilità di una condizione.

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

Avviare ShinyStan per la diagnostica interattiva

shinystan::launch_shinystan(fit) apre un'app Shiny interattiva con traceplot, distribuzioni a posteriori, grafici pairs e diagnostica NUTS, tutto in un unico ambiente. È lo strumento più completo per esplorare un fit 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

Correzioni comuni per la convergenza

Quando Rhat > 1.01 o si osservano divergenze:

Aumenti iter e warmup
Aumenti adapt_delta verso 1.0 (ad esempio, 0.95) in control
Riparametrizzi il modello — usi una parametrizzazione non centrata per i modelli gerarchici
Restringa le prior se sono troppo diffuse
Controlli i dati — i valori anomali o le differenze di scala causano problemi al sampler

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

Verifica rapida: soglia Rhat

Qual è la soglia standard attuale di Rhat che indica che un modello Stan ha raggiunto la convergenza?

Riepilogo del campionamento e della diagnostica MCMC

Flusso di lavoro MCMC in RStan:

stan(model_code=..., data=..., chains=4, iter=2000, warmup=1000) adatta il modello
print(fit) mostra Rhat e n_eff — i principali diagnostici di convergenza
Rhat < 1.01 e n_eff / total > 0.1 indicano un campione con comportamento adeguato
traceplot() — verifica visiva del mixing; pairs() — rivela problemi nella geometria a posteriori
extract(fit, pars='mu')$mu — consente di accedere ai campioni a posteriori grezzi
shinystan::launch_shinystan(fit) — diagnostica interattiva completa

Domande Frequenti

La lezione «Campionamento MCMC e diagnostica» è gratuita?

Sì — il testo completo di «Campionamento MCMC e diagnostica» è gratuito qui sul web. Per esercitarvi in modo interattivo (un editor di codice integrato e un tutor IA 24/7) e sbloccare il resto del corso R Academy, passa a CoddyKit PRO. Il corso R Academy include 4 lezioni in totale.

Cosa imparerò in «Campionamento MCMC e diagnostica»?

Esegua il campionamento, esamini le catene e interpreti le diagnostiche Rhat ed ESS Eserciti R Academy con codice pratico che esegui direttamente nel browser, e un tutor IA 24/7 risponde alle tue domande mentre lavori sulla lezione.

Ho bisogno di esperienza per iniziare R Academy?

Non è richiesta alcuna esperienza precedente. R Academy su CoddyKit è strutturato per principianti e studenti avanzati, quindi puoi iniziare da qui o dall'inizio e procedere al tuo ritmo. Questa è la lezione 3 di 4.

Quanto tempo richiede la lezione «Campionamento MCMC e diagnostica»?

La maggior parte delle lezioni CoddyKit richiede circa 5–10 minuti. Ogni lezione è breve e interattiva, quindi fai progressi costanti e riprendi esattamente da dove hai lasciato su web e app.

Posso scrivere ed eseguire codice in questa lezione R Academy?

Sì. Ogni lezione R Academy include un editor di codice integrato, quindi scrivi ed esegui codice reale direttamente nel tuo browser e ricevi feedback istantaneo dall'IA — nessuna configurazione locale necessaria.

Tutte le lezioni di questo corso

← Torna a R Academy