here::here("code", "_common.R") |>
source()
# Load packages
if (!requireNamespace("pacman")) install.packages("pacman")
pacman::p_load(cmdstanr, reshape2)
56 L’algoritmo di Metropolis-Hastings
Introduzione
Nei capitoli precedenti abbiamo esaminato diversi esempi di inferenza bayesiana, concentrandoci su modelli parametrici semplici come il caso unidimensionale del modello bernoulliano. In questi scenari elementari, abbiamo utilizzato due approcci principali: l’approssimazione discreta su griglia, che valuta la distribuzione a posteriori su un insieme finito di valori parametrali, e l’uso di distribuzioni a priori coniugate, che, grazie alla loro forma analiticamente trattabile, permettono di derivare la posterior in modo esatto.
Tuttavia, questi metodi presentano limitazioni significative non appena ci si allontana da contesti idealizzati. In particolare, diventano rapidamente inefficienti in presenza di modelli con molteplici parametri, relazioni complesse tra variabili, o distribuzioni a priori non coniugate. In tali casi, il calcolo analitico risulta spesso intrattabile, mentre l’approssimazione su griglia diventa computazionalmente proibitiva a causa della crescita esponenziale del numero di valutazioni necessarie al crescere della dimensionalità dello spazio parametrico.
Per affrontare queste sfide, l’inferenza bayesiana moderna fa ricorso a tecniche numeriche avanzate, tra cui spicca il Metodo Monte Carlo a Catena di Markov (MCMC). Questo approccio innovativo supera i limiti dei metodi tradizionali generando campioni dalla distribuzione a posteriori attraverso un processo iterativo che combina l’efficienza del campionamento stocastico con la flessibilità degli algoritmi markoviani. Il risultato è una metodologia potente e generale, in grado di fornire stime accurate anche per modelli complessi e ad alta dimensionalità, senza richiedere soluzioni analitiche chiuse.
56.1 L’obiettivo del metodo MCMC
Il metodo MCMC è un approccio computazionale che consente di approssimare distribuzioni di probabilità complesse, generando una sequenza di valori campionati che seguono la distribuzione a posteriori di interesse.
L’idea di base è la seguente:
- consideriamo la distribuzione a posteriori come una popolazione da cui desideriamo estrarre campioni;
- generando un numero sufficientemente grande di campioni (ad esempio, diverse migliaia), la distribuzione empirica dei campioni ottenuti si avvicina progressivamente alla distribuzione teorica a posteriori;
- in questo modo, possiamo stimare quantità di interesse, come la media, la varianza, o intervalli di credibilità, anche senza conoscere una forma analitica esplicita della distribuzione a posteriori.
56.1.1 La natura dipendente del campionamento MCMC
A differenza delle tecniche di campionamento indipendente precedentemente esaminate, l’approccio MCMC genera una sequenza di valori correlati attraverso un meccanismo di transizione markoviana. La caratteristica distintiva di questo processo risiede nella proprietà di Markov: ogni nuovo campione dipende esclusivamente dallo stato corrente della catena, mostrando memoria soltanto a breve termine piuttosto che dipendenza dall’intera storia precedente.
Questa architettura sequenziale produce inevitabilmente autocorrelazione tra osservazioni adiacenti. Quando la catena visita una regione ad alta densità della distribuzione target, tenderà a persistere in tale zona per diverse iterazioni prima di migrare verso altre regioni. Mentre tale comportamento è funzionale all’esplorazione efficiente dello spazio parametrico, introduce importanti considerazioni pratiche:
- l’informazione effettiva contenuta in
campioni correlati è inferiore a quella di campioni indipendenti; - la valutazione della convergenza richiede analisi diagnostiche specifiche;
- la dimensione efficace del campione (ESS) diventa un parametro cruciale per la qualità dell’inferenza.
Per compensare questa riduzione di informazione per campione, è generalmente necessario generare sequenze più lunghe rispetto al campionamento indipendente. Tuttavia, questo svantaggio apparente è ampiamente compensato dalla capacità di MCMC di affrontare problemi complessi che risulterebbero intrattabili con metodi tradizionali.
56.1.2 Perché utilizzare MCMC
Il metodo MCMC è diventato uno strumento centrale nell’inferenza bayesiana contemporanea perché:
- è in grado di affrontare problemi complessi, caratterizzati da distribuzioni a posteriori di forma irregolare o definite in spazi ad alta dimensione;
- non richiede il calcolo diretto dell’integrale di normalizzazione che compare nel teorema di Bayes;
- permette di ottenere approssimazioni accurate della distribuzione a posteriori tramite simulazione numerica.
Nel seguito ci concentreremo sull’algoritmo di Metropolis, uno dei metodi più semplici ed essenziali per implementare il campionamento MCMC.
56.2 L’algoritmo di Metropolis: introduzione intuitiva
L’algoritmo di Metropolis è un metodo MCMC che consente di esplorare una distribuzione di probabilità complessa costruendo una sequenza di campioni dipendenti tra loro.
La logica dell’algoritmo può essere riassunta nei seguenti passaggi fondamentali:
-
Punto di partenza: si inizia da un valore iniziale
scelto arbitrariamente. -
Proposta di un nuovo punto: si genera un nuovo valore candidato
partendo da , utilizzando una distribuzione di proposta (ad esempio una distribuzione normale centrata su ). -
Valutazione della proposta: si confrontano le densità a posteriori associate al valore attuale
e al valore proposto . -
Decisione di accettazione:
- se
ha una densità a posteriori più alta di , viene accettato automaticamente; - se
ha una densità a posteriori inferiore, viene accettato con una certa probabilità proporzionale al rapporto delle densità.
- se
- Registrazione: in ogni caso, si registra la posizione attuale (sia che si sia accettato un nuovo punto, sia che si sia rimasti fermi).
Questo processo viene ripetuto per un numero elevato di iterazioni, generando una catena di campioni che, dopo un opportuno periodo iniziale (detto burn-in), approssima la distribuzione a posteriori.
56.3 Perché accettiamo anche mosse peggiori
Uno degli aspetti peculiari dell’algoritmo di Metropolis è la possibilità di accettare anche proposte peggiori, ossia punti
Questa scelta ha una motivazione fondamentale:
- se accettassimo solo le mosse che migliorano la densità, l’algoritmo rischierebbe di bloccarsi in un massimo locale della distribuzione, senza esplorare altre aree che, pur avendo densità più bassa localmente, potrebbero condurre a regioni più interessanti globalmente;
- accettare occasionalmente mosse peggiori consente all’algoritmo di esplorare meglio tutto lo spazio dei parametri, evitando di rimanere intrappolato in una singola area.
In questo modo, la catena può attraversare regioni di bassa probabilità e raggiungere altre modalità della distribuzione, garantendo una copertura più completa dello spazio delle soluzioni plausibili.
56.4 La scelta della larghezza della proposta
Nell’algoritmo di Metropolis, la proposta di un nuovo valore
La scelta del valore di
- se
è troppo piccolo, i passi proposti saranno molto vicini al punto attuale. In questo caso, molte proposte saranno accettate, ma la catena si muoverà lentamente nello spazio dei parametri, esplorandolo inefficientemente (alta correlazione tra i campioni); - se
è troppo grande, i passi proposti saranno molto lontani dal punto attuale. In questo caso, la maggior parte delle proposte cadrà in regioni di bassa densità, portando a un alto tasso di rifiuto delle proposte e quindi a una scarsa efficienza del campionamento.
Un valore ottimale di
- accettazione sufficiente di nuove proposte;
- esplorazione efficiente dello spazio dei parametri.
In generale, si cerca di ottenere un tasso di accettazione compreso tra il 40% e il 50% per l’algoritmo di Metropolis a singolo parametro.
56.5 L’importanza dei grafici diagnostici: Trace plot e Correlogramma
Per valutare la qualità della catena generata dall’algoritmo di Metropolis, è fondamentale analizzare alcuni grafici diagnostici.
56.5.1 Trace plot
Il trace plot rappresenta i valori campionati di
Un trace plot di buona qualità mostra:
- oscillazioni attorno a un valore centrale stabile (assenza di trend sistematici);
- una copertura adeguata dello spazio plausibile per
.
Un trace plot problematico può rivelare:
- fasi iniziali instabili (burn-in non sufficientemente lungo);
- mancata esplorazione completa della distribuzione;
- convergenza solo apparente, con la catena bloccata in una modalità.
56.5.2 Correlogramma
Il correlogramma mostra il grado di autocorrelazione dei campioni in funzione del numero di passi di lag.
Idealmente:
- l’autocorrelazione dovrebbe decrescere rapidamente all’aumentare del lag;
- una catena ben mescolata presenta un correlogramma che si avvicina rapidamente a zero.
Una forte autocorrelazione indica che i campioni successivi sono troppo simili tra loro, riducendo l’efficienza dell’inferenza statistica.
Questi concetti costituiscono il fondamento necessario per affrontare la comprensione operativa e pratica dell’algoritmo di Metropolis che svilupperemo nei prossimi esempi. A questo fine, il capitolo è strutturato in varie sezioni che facilitano la comprensione progressiva del tema.
- Inizieremo discutendo di come la distribuzione a posteriori possa essere approssimata mediante tecniche di simulazione convenzionali. Questa prima parte presuppone che la distribuzione target, o “a posteriori,” sia già conosciuta o disponibile per l’analisi.
- In seguito, passeremo a illustrare come l’algoritmo di Metropolis possa essere utilizzato per affrontare situazioni in cui la distribuzione a posteriori non è direttamente nota. In questi casi, spesso abbiamo a disposizione informazioni riguardanti la distribuzione a priori e la funzione di verosimiglianza, che possono essere utilizzate per ottenere un’approssimazione efficace della distribuzione a posteriori.
56.6 Un Esempio Concreto
A titolo esemplificativo, utilizzeremo il dataset moma_sample.csv
, il quale costituisce un campione casuale di 100 artisti provenienti dal Museo di Arte Moderna di New York (MoMA) e contiene diverse informazioni relative a ciascun artista.
Il nostro interesse è focalizzato sulla determinazione della probabilità che un artista presente nel MoMA appartenga alla generazione X o a una generazione successiva (nati dopo il 1965). Questa probabilità sarà indicata come
Importiamo i dati.
Esaminiamo le prime cinque righe del DataFrame.
moma_sample |>
head()
#> artist country birth death alive genx gender count
#> 1 Ad Gerritsen dutch 1940 2015 FALSE FALSE male 1
#> 2 Kirstine Roepstorff danish 1972 NA TRUE TRUE female 3
#> 3 Lisa Baumgardner american 1958 2015 FALSE FALSE female 2
#> 4 David Bates american 1952 NA TRUE FALSE male 1
#> 5 Simon Levy american 1946 NA TRUE FALSE male 1
#> 6 Pierre Mercure canadian 1927 1966 FALSE FALSE male 8
#> year_acquired_min year_acquired_max
#> 1 1981 1981
#> 2 2005 2005
#> 3 2016 2016
#> 4 2001 2001
#> 5 2012 2012
#> 6 2008 2008
Dai dati osserviamo che solo 14 artisti su 100 appartengono alla generazione X o a una generazione successiva.
# Calcoliamo la distribuzione delle generazioni
result <- table(moma_sample$genx)
result
#>
#> FALSE TRUE
#> 86 14
Il valore campionato
Possiamo interpretare i dati
Supponiamo che le nostre credenze pregresse riguardo a
Sfruttando le proprietà delle distribuzioni coniugate, possiamo calcolare esattamente la distribuzione a posteriori:
# Y ~ Binomiale(100, π)
# θ ~ Beta(4, 6)
# Posteriori: θ | (Y = 14) ~ Beta(4 + 14, 6 + 100 - 14) → Beta(18, 92)
Nella figura seguente, è rappresentata la distribuzione a posteriori del parametro
# Sequenza per θ
x <- seq(0, 1, length.out = 1000)
# Densità prior e posterior
prior_density <- dbeta(x, 4, 6)
posterior_density <- dbeta(x, 18, 92)
# Dati per il grafico
df <- data.frame(
x = x,
prior = prior_density,
posterior = posterior_density
)
# Lungo
df_long <- reshape2::melt(
df,
id.vars = "x",
measure.vars = c("prior", "posterior"),
variable.name = "distribuzione",
value.name = "densita"
)
# Fattorizziamo con levels e labels espliciti
df_long <- dplyr::mutate(
df_long,
distribuzione = factor(
distribuzione,
levels = c("prior", "posterior"),
labels = c("Prior: Beta(4, 6)", "Posterior: Beta(18, 92)")
)
)
if (!exists("palette_set1")) {
palette_set1 <- RColorBrewer::brewer.pal(3, "Set1")
names(palette_set1) <- c("uno","due","tre")
}
cols <- c(
"Prior: Beta(4, 6)" = unname(palette_set1[[1]]),
"Posterior: Beta(18, 92)" = unname(palette_set1[[2]])
)
# Grafico
ggplot(df_long, aes(x = x, y = densita)) +
geom_line(aes(color = distribuzione), linewidth = 1) +
geom_area(aes(fill = distribuzione), alpha = 0.35, position = "identity") +
scale_color_manual(values = cols, breaks = names(cols), labels = names(cols)) +
scale_fill_manual(values = cols, breaks = names(cols), labels = names(cols)) +
labs(
title = "Densità a priori e a posteriori",
x = "Valore del parametro",
y = "Densità",
color = NULL, fill = NULL
) +
theme(
plot.title = element_text(hjust = 0.5),
legend.position = "top"
)
In questo grafico, la curva blu rappresenta la distribuzione a priori
Se vogliamo conoscere il valore della media a posteriori di
56.6.1 Simulazione con distribuzione target nota
Usiamo ora una simulazione numerica per stimare la media a posteriori di
Se vogliamo ottenere un risultato approssimato con un numero limitato di campioni (ad esempio, 10), possiamo utilizzare la seguente simulazione:
# Calcoliamo la media dei campioni
mean(y)
#> [1] 0.1508
Tuttavia, con soli 10 campioni, l’approssimazione potrebbe non essere molto accurata. Aumentando il numero di campioni, ad esempio a 10,000, possiamo ottenere una stima molto più precisa:
Quando il numero di campioni dalla distribuzione a posteriori diventa molto grande, la media campionaria converge al valore atteso della distribuzione della popolazione. Questo principio non si applica solo alla media, ma anche ad altre statistiche descrittive come la moda e la varianza.
È importante sottolineare che l’applicazione della simulazione di Monte Carlo è efficace per calcolare distribuzioni a posteriori solo quando conosciamo la distribuzione stessa e possiamo utilizzare funzioni R per estrarre campioni casuali da tale distribuzione. Ciò è stato possibile nel caso della distribuzione a posteriori
Tuttavia, questa situazione ideale non si verifica sempre nella pratica, poiché le distribuzioni a priori coniugate alla verosimiglianza sono spesso rare. Per esempio, nel caso di una verosimiglianza binomiale e una distribuzione a priori gaussiana, l’espressione
rende impossibile calcolare analiticamente la distribuzione a posteriori di
In queste circostanze, però, è possibile ottenere campioni casuali dalla distribuzione a posteriori mediante l’uso di metodi Monte Carlo basati su Catena di Markov (MCMC). Gli algoritmi MCMC, come ad esempio l’algoritmo Metropolis, costituiscono una classe di metodi che consentono di estrarre campioni casuali dalla distribuzione a posteriori senza richiedere la conoscenza della sua rappresentazione analitica. Le tecniche MCMC sono ampiamente adottate per risolvere problemi di inferenza bayesiana e rappresentano il principale strumento computazionale per ottenere stime approssimate di distribuzioni a posteriori in situazioni complesse e non analiticamente trattabili.
56.6.2 Algoritmo di Metropolis
L’algoritmo di Metropolis appartiene alla famiglia dei metodi Monte Carlo basati su catene di Markov (MCMC), sfruttando le proprietà di queste catene per generare campioni da una distribuzione target. Il suo obiettivo principale è di esplorare lo spazio dei parametri in modo efficiente, producendo campioni che approssimano la distribuzione a posteriori di interesse.
56.6.3 Principio di Funzionamento
L’algoritmo inizia da un valore iniziale per i parametri e, in ogni iterazione, genera un nuovo campione tramite una distribuzione di proposta (solitamente una distribuzione normale centrata sul valore corrente). Successivamente, decide se accettare il nuovo campione in base al confronto tra le densità posteriori del nuovo campione e di quello precedente. Questa accettazione avviene in modo probabilistico, favorendo campioni con una densità più alta ma consentendo anche l’accettazione di campioni peggiori per evitare che la catena rimanga bloccata in minimi locali.
56.6.4 Burn-in e Convergenza
Poiché i primi campioni potrebbero non rappresentare bene la distribuzione target, si esclude spesso una porzione iniziale della catena (fase di burn-in). Con il progredire delle iterazioni, i campioni si distribuiscono in accordo con la distribuzione stazionaria desiderata, indipendentemente dallo stato iniziale scelto. Questo processo permette di esplorare lo spazio dei parametri in modo efficiente.
56.6.5 Meccanismo di Accettazione e Rifiuto
L’algoritmo di Metropolis bilancia due esigenze opposte:
- Esplorazione di nuove aree dello spazio dei parametri.
- Sfruttamento delle informazioni già acquisite dai campioni precedenti.
Utilizzando una regola probabilistica per accettare campioni peggiori (con minore densità a posteriori), l’algoritmo evita di restare intrappolato in minimi locali, esplorando così in modo più completo l’intera distribuzione.
56.6.6 Passaggi Fondamentali dell’Algoritmo di Metropolis
-
Scelta di uno stato iniziale
e impostazione del contatore .- Questo è il punto di partenza della catena, dove
rappresenta il primo campione.
- Questo è il punto di partenza della catena, dove
-
Proposta di un nuovo campione
.- Un nuovo valore
viene generato da una distribuzione di proposta , solitamente una distribuzione normale centrata sul campione corrente con una deviazione standard che controlla l’ampiezza dei passi.
- Un nuovo valore
-
Verifica dei vincoli del campione proposto.
- Se il nuovo campione deve rispettare dei vincoli (ad esempio, essere compreso tra 0 e 1 per probabilità), campioni non validi vengono automaticamente rifiutati.
-
Calcolo del rapporto di accettazione
.- Si calcola
, che rappresenta il rapporto tra le densità a posteriori del nuovo campione e del campione corrente . Questo valore guida la decisione di accettazione.
- Si calcola
-
Decisione di accettazione.
- Se
, il nuovo campione viene accettato incondizionatamente. - Se
, il campione viene accettato con probabilità . In caso di rifiuto, si mantiene il campione corrente per la prossima iterazione.
- Se
-
Ripetizione del processo.
- Si ripetono i passaggi dal 2 al 5 fino a ottenere il numero desiderato di campioni.
56.6.7 Dettagli Aggiuntivi
Distribuzione di proposta: La distribuzione di proposta
genera nuovi campioni attorno a . Tipicamente si usa una normale , dove controlla quanto il nuovo campione si discosta da quello corrente. Scegliere un troppo piccolo può rendere l’esplorazione lenta, mentre un troppo grande può far rifiutare troppi campioni, riducendo l’efficienza.Rapporto di accettazione
: Se il nuovo campione ha una densità a posteriori maggiore del campione corrente, viene sempre accettato. Se ha una densità inferiore, viene accettato con probabilità , il che consente di esplorare anche regioni meno probabili della distribuzione.Accettazione probabilistica: Accettare campioni peggiori occasionalmente aiuta l’algoritmo a evitare di bloccarsi in minimi locali. Questo è uno dei punti di forza dell’algoritmo di Metropolis, che garantisce una buona esplorazione dello spazio dei parametri.
56.7 Esempio di Implementazione
Supponiamo di voler stimare la probabilità
56.7.1 Definizione della Distribuzione a Priori
La funzione prior
calcola la densità della distribuzione Beta(4, 6) per un dato
# Distribuzione a priori Beta(4, 6)
prior <- function(p) {
dbeta(p, shape1 = 4, shape2 = 6)
}
Questa distribuzione esprime la nostra plausibilità iniziale sui valori di
56.7.2 Funzione di Verosimiglianza
La funzione likelihood
modella la probabilità di osservare 14 successi su 100 prove:
# Verosimiglianza binomiale (y = 14 successi su n = 100 prove)
likelihood <- function(p) {
y <- 14
n <- 100
dbinom(y, size = n, prob = p)
}
56.7.3 Distribuzione a Posteriori Non Normalizzata
La posteriori si ottiene combinando priori e verosimiglianza:
# Posteriori non normalizzata (prodotto tra verosimiglianza e priori)
posterior <- function(p) {
likelihood(p) * prior(p)
}
56.7.4 Distribuzione Proposta
La distribuzione proposta sarà una distribuzione normale centrata sullo stato corrente con una deviazione standard specificata:
# Generazione proposta (normale con media sullo stato corrente)
proposal_distribution <- function(current_state, proposal_sigma) {
rnorm(1, mean = current_state, sd = proposal_sigma)
}
56.7.5 Implementazione dell’Algoritmo Metropolis-Hastings
Procediamo ora con l’implementazione dell’algoritmo di Metropolis-Hastings, considerando i dati relativi agli artisti della Generazione X presenti al MoMA. La distribuzione a priori per
# Algoritmo Metropolis-Hastings
metropolis_hastings <- function(n_samples, start, proposal_sigma) {
samples <- numeric(n_samples)
current <- start # Stato iniziale
for (i in seq_len(n_samples)) {
proposal <- proposal_distribution(current, proposal_sigma)
# Verifica validità e calcolo rapporto di accettazione
if (proposal >= 0 && proposal <= 1) {
acceptance_ratio <- min(1, posterior(proposal) / posterior(current))
# Accetta/rifiuta con probabilità acceptance_ratio
if (runif(1) < acceptance_ratio) {
current <- proposal
}
}
samples[i] <- current # Aggiorna la catena
}
return(samples)
}
56.7.6 Esecuzione dell’Algoritmo
# Parametri dell'algoritmo
n_samples <- 10000
start <- 0.5
proposal_sigma <- 0.1
# Esecuzione del campionamento
set.seed(123) # Per riproducibilità
samples <- metropolis_hastings(n_samples, start, proposal_sigma)
56.7.7 Analisi dei Risultati
Scartiamo i primi 5000 campioni per considerare solo quelli generati dopo il burn-in:
Calcoliamo la media e la deviazione standard dei campioni:
Visualizziamo l’evoluzione della catena per i primi 200 campioni e per quelli post-burn-in:
tibble(
Iterazione = 1:200,
Theta = samples[1:200]
) |>
ggplot(aes(x = Iterazione, y = Theta)) +
geom_line() + # Specifica che vuoi un grafico a linea
ggtitle("Trace Plot (Primi 200 Campioni)") + # Aggiunge il titolo principale
xlab("Iterazioni") + # Etichetta l'asse X
ylab(expression(theta)) # Etichetta l'asse Y usando l'espressione per theta
tibble(
Iterazione = 1:length(post_burnin_samples),
Theta = post_burnin_samples
) |>
ggplot(aes(x = Iterazione, y = Theta)) +
geom_line() + # Specifica un grafico a linea
ggtitle("Trace Plot (Post Burn-in)") + # Aggiunge il titolo
xlab("Iterazioni") + # Indice all'interno della serie post-burn-in
ylab(expression(theta)) # Etichetta l'asse Y
Sovrapponiamo la distribuzione analitica
# Se per qualche motivo non è stato caricato _common.R:
if (!exists("palette_set1")) {
palette_set1 <- RColorBrewer::brewer.pal(3, "Set1")
names(palette_set1) <- c("uno", "due", "tre")
}
# Mappa colori coerente con la figura precedente
colori_grafico <- c(
"Istogramma MCMC" = unname(palette_set1[[1]]),
"Beta(18, 92)" = unname(palette_set1[[2]])
)
tibble(Theta = post_burnin_samples) |>
ggplot(aes(x = Theta)) +
geom_histogram(
aes(y = after_stat(density), fill = "Istogramma MCMC"),
bins = 30,
color = "black",
alpha = 0.7
) +
stat_function(
aes(color = "Beta(18, 92)"),
fun = dbeta,
args = list(shape1 = 18, shape2 = 92),
linewidth = 1.2
) +
labs(
title = "Istogramma e distribuzione a posteriori",
x = expression(theta),
y = "Densità",
fill = "Distribuzione",
color = "Distribuzione"
) +
scale_fill_manual(values = colori_grafico, breaks = names(colori_grafico)) +
scale_color_manual(values = colori_grafico, breaks = names(colori_grafico)) +
theme(
plot.title = element_text(hjust = 0.5),
legend.position = "top"
)
Calcoliamo l’intervallo di credibilità al 94%:
I valori ottenuti con l’algoritmo di Metropolis (usando solo un piccolo numero di iterazioni) sono quasi identici ai valori esatti:
Questa implementazione in R dimostra come utilizzare l’algoritmo di Metropolis per stimare una distribuzione a posteriori e analizzare i risultati in modo dettagliato e riproducibile.
56.8 Catene di Markov e Convergenza
Nell’ambito delle simulazioni Monte Carlo, una catena rappresenta una sequenza di valori campionati dall’algoritmo durante le sue iterazioni. Ogni valore nella catena corrisponde a un possibile stato del sistema che stiamo modellando. In altre parole, una catena traccia il percorso che l’algoritmo segue nello spazio dei parametri, esplorando le diverse configurazioni possibili.
Per verificare se l’algoritmo ha raggiunto la convergenza e se i campioni generati rappresentano effettivamente la distribuzione di interesse, è utile eseguire multiple catene. Ogni catena parte da un punto iniziale diverso nello spazio dei parametri.
I vantaggi delle multiple catene:
- Diagnostica della convergenza: Confrontando le diverse catene, possiamo valutare se si stabilizzano verso la stessa distribuzione. Se le catene si mescolano bene, ovvero si intersecano frequentemente nel grafico dei valori campionati (trace plot), è un forte indicatore di convergenza.
- Robustezza: L’utilizzo di multiple catene rende l’analisi meno sensibile alla scelta del punto di partenza. Se una singola catena potesse rimanere “intrappolata” in una regione dello spazio dei parametri, multiple catene aumentano la probabilità di esplorare lo spazio in modo più completo.
56.9 Diagnostiche della soluzione MCMC
56.9.1 Stazionarietà e Convergenza
Un aspetto cruciale nell’analisi delle catene di Markov MCMC è la convergenza alla distribuzione stazionaria. Intuitivamente, la catena converge quando i campioni generati rappresentano fedelmente la distribuzione di interesse, indipendentemente dal punto di partenza. Questo fenomeno è spesso indicato come “mixing”.
56.9.1.1 Valutazione Visuale: Trace Plots e Grafici di Densità
- Trace Plots: Questi grafici visualizzano l’evoluzione dei parametri nel tempo. Una catena convergente mostra tracce stabili e senza trend evidenti. Tracce irregolari o con andamenti sistematici suggeriscono problemi di convergenza.
- Grafici di Densità: Confrontando i grafici di densità dei campioni con la distribuzione teorica, è possibile valutare visivamente se la catena sta esplorando adeguatamente lo spazio dei parametri. Una buona convergenza si manifesta con una sovrapposizione tra i due grafici.
Segni di Convergenza:
- Stabilità: I valori campionati oscillano attorno a un valore medio costante, senza trend marcati.
- Omogeneità: La variabilità dei campioni rimane relativamente uniforme nel tempo.
- Assenza di Periodicità: Non si osservano pattern ciclici o ripetitivi.
In sintesi, i trace plots e i grafici di densità offrono strumenti visivi rapidi per valutare la convergenza di una catena di Markov MCMC. Una convergenza soddisfacente è fondamentale per garantire la validità delle inferenze statistiche basate sui campioni generati.
56.9.2 Autocorrelazione nelle catene di Markov MCMC
A differenza dei generatori di numeri casuali indipendenti, gli algoritmi MCMC producono una sequenza di campioni correlati. Ogni valore campionato dipende da quello precedente, formando una catena di Markov. Questa interdipendenza è un aspetto fondamentale dell’MCMC.
L’autocorrelazione quantifica il grado di dipendenza tra valori distanti di una certa quantità (detta lag) nella catena. Un’alta autocorrelazione a lag bassi indica una forte dipendenza tra campioni successivi. Al contrario, una rapida diminuzione dell’autocorrelazione al crescere del lag suggerisce che la catena “miscela” bene, ovvero esplora lo spazio dei parametri in modo efficiente.
- Lag 1: Misura la correlazione tra valori consecutivi nella catena.
- Lag 2: Misura la correlazione tra valori separati da un passo intermedio.
- Lag k: Generalizza il concetto ai valori separati da k passi.
Un correlogramma è un grafico che mostra l’autocorrelazione in funzione del lag. Un decadimento rapido dell’autocorrelazione verso zero indica una buona convergenza della catena.
L’autocorrelazione di ordine
56.9.3 Esempio di Simulazione di Dati Autocorrelati
Per fare un esempio pratico, creiamo un vettore di dati autocorrelati:
# Creiamo un vettore di dati
x <- c(22, 24, 25, 25, 28, 29, 34, 37, 40, 44, 51, 48, 47, 50, 51)
x
#> [1] 22 24 25 25 28 29 34 37 40 44 51 48 47 50 51
56.9.3.1 Calcolo dell’Autocorrelazione
L’autocorrelazione di ordine 1 è la correlazione tra ciascun elemento e il successivo nella sequenza. In R possiamo utilizzare la funzione acf()
per calcolare l’autocorrelazione.
# Calcolo dell'autocorrelazione
acf_values <- acf(x, plot = FALSE)
acf_values
#>
#> Autocorrelations of series 'x', by lag
#>
#> 0 1 2 3 4 5 6 7 8 9
#> 1.000 0.832 0.656 0.491 0.279 0.031 -0.165 -0.304 -0.401 -0.458
#> 10 11
#> -0.450 -0.369
Nell’esempio, il vettore x
rappresenta una serie temporale di 15 elementi. Il calcolo dell’autocorrelazione restituisce i seguenti valori per i primi ritardi (lag):
- 0.8317: autocorrelazione di ordine 1 (lag = 1),
- 0.6563: autocorrelazione di ordine 2 (lag = 2),
-
0.4910: autocorrelazione di ordine 3 (lag = 3),
ecc.
56.9.3.2 Specifica del Numero di Ritardi (Lag)
Possiamo limitare il numero di ritardi calcolati utilizzando l’argomento lag.max
nella funzione acf()
:
# Calcolo dell'autocorrelazione per i primi 4 lag
acf(x, lag.max = 4, plot = FALSE)
#>
#> Autocorrelations of series 'x', by lag
#>
#> 0 1 2 3 4
#> 1.000 0.832 0.656 0.491 0.279
56.9.3.3 Grafico della Funzione di Autocorrelazione (Correlogramma)
In R possiamo creare un correlogramma con la funzione acf()
:
# Correlogramma per la serie temporale
acf(x, main = "Correlogramma della Serie Temporale", lag.max = 9)
56.9.4 Analisi della Catena di Markov
Applichiamo lo stesso approccio alla catena di Markov ottenuta precedentemente, considerando i campioni post burn-in:
In situazioni ideali, l’autocorrelazione diminuisce rapidamente, diventando insignificante per piccoli lag. Questo comportamento è un’indicazione del “mixing” efficace della catena, ossia della sua convergenza alla distribuzione stazionaria.
56.9.5 Sottocampionamento (Thinning)
Per ridurre l’autocorrelazione, possiamo applicare una strategia di sottocampionamento (thinning), memorizzando solo ogni
In conclusione, il correlogramma con thinning mostra che l’autocorrelazione diminuisce più rapidamente rispetto ai campioni originali, suggerendo che la strategia di sottocampionamento è efficace nel migliorare l’indipendenza tra i campioni successivi. Questo migliora la qualità delle inferenze basate sulla catena di Markov.
56.9.5.1 Tasso di accettazione
Quando si utilizza l’algoritmo Metropolis, è importante monitorare il tasso di accettazione e assicurarsi che sia nell’intervallo ottimale. Se si accetta quasi sempre il candidato proposto, probabilmente significa che, in ogni iterazione, la catena salta solo di un piccolo passo (in modo che il rapporto di accettazione sia vicino a 1 ogni volta). Di conseguenza, la catena impiegherà molte iterazioni per raggiungere altre regioni della distribuzione stazionaria e i campioni consecutivi saranno molto fortemente correlati. D’altra parte, se il tasso di accettazione è molto basso, la catena rimarrà bloccata nella stessa posizione per molte iterazioni prima di spostarsi verso uno stato diverso. Per l’algoritmo Metropolis base con un singolo parametro con una distribuzione proposta Gaussiana normale, un tasso di accettazione ottimale è compreso tra il 40% e il 50%.
56.9.6 Test Statistici per la Convergenza
Oltre agli approcci grafici, esistono test statistici specifici che possono aiutare a determinare se la catena ha raggiunto uno stato stazionario.
56.9.6.1 Test di Geweke
Il test di Geweke è una procedura che confronta le medie di due segmenti della catena di campionamento, tipicamente il primo 10% e l’ultimo 50% dei campioni, dopo aver escluso un iniziale periodo di “burn-in” (una fase iniziale durante la quale la catena potrebbe non essere ancora convergente). La premessa di base è che, se la catena è in uno stato stazionario, le medie di questi due segmenti dovrebbero essere sostanzialmente uguali. Differenze importanti tra queste medie possono indicare che la catena non ha ancora raggiunto la convergenza.
56.9.6.2 Geweke Z-score
Una variante del test di Geweke è lo z-score di Geweke, che offre un modo quantitativo per valutare le differenze tra i segmenti della catena. Questo test calcola uno z-score che confronta le medie dei due segmenti tenendo conto della varianza. Un valore di z-score:
- Al di sotto di 2 (in valore assoluto) suggerisce che non ci sono differenze degne di nota tra i segmenti, indicando che la catena potrebbe essere in stato stazionario.
- Superiore a 2 (in valore assoluto) indica che esiste una differenza degna di nota tra i segmenti, suggerendo che la catena non ha raggiunto la convergenza e potrebbe essere necessario un periodo di burn-in più esteso.
Entrambi i metodi forniscono strumenti utili per valutare la convergenza delle catene MCMC. È importante notare che nessun test può garantire con certezza la convergenza, ma l’utilizzo congiunto di approcci grafici e test statistici può offrire una buona indicazione dello stato della catena.
56.9.7 Dimensione del campione effettiva (ESS)
La correlazione tra campioni consecutivi in una catena MCMC riduce l’informazione effettiva contenuta in ogni iterazione. La dimensione del campione effettiva (ESS) quantifica questa perdita di informazione dovuta alla dipendenza tra i campioni, stimando il numero equivalente di campioni indipendenti. Un valore basso di ESS indica una forte correlazione tra i campioni e una convergenza più lenta della catena.
L’ESS descrive l’efficacia del campionamento dipendente in termini di campioni indipendenti estratti dalla stessa distribuzione. Rappresenta un indicatore dell’efficienza del campionamento e dell’autocorrelazione della catena.
La formula per stimare la dimensione del campione effettiva (ESS) di una catena di Markov è:
dove:
-
è il numero totale di campioni nella catena, -
è il lag, ovvero il numero massimo di termini di autocorrelazione considerati, -
è l’autocorrelazione al lag , ossia la correlazione tra due campioni consecutivi separati da iterazioni.
In pratica,
56.9.8 Calcolo della Statistica di Gelman-Rubin ( )
Per calcolare la statistica di Gelman-Rubin (spesso indicata come
- Esegui
catene di Markov di lunghezza , dove è solitamente maggiore di 1. - Per ciascun parametro scalare
, calcola la varianza all’interno delle catene ( ) e la varianza tra le catene ( ). - Calcola la varianza combinata
come media ponderata delle varianze all’interno delle catene. - Calcola il fattore di riduzione della scala potenziale
come la radice quadrata del rapporto tra la varianza combinata e la varianza all’interno delle catene :
- Se
è vicino a 1, ciò indica che le catene sono in convergenza.
La statistica di Gelman-Rubin
56.10 Vantaggi del Campionamento MCMC rispetto alle Soluzioni Analitiche
Il campionamento MCMC offre notevoli vantaggi pratici rispetto alle soluzioni analitiche nella statistica bayesiana, in particolare quando si tratta di manipolare distribuzioni a posteriori. Sebbene l’impossibilità di calcolare analiticamente la distribuzione a posteriori sia spesso la motivazione principale per l’uso di MCMC, i benefici di questo approccio si estendono ben oltre questa necessità (Bürkner, 2024).
56.10.1 Facilità di Manipolazione e Flessibilità
Il vantaggio chiave del campionamento MCMC risiede nella semplicità con cui si possono manipolare i campioni ottenuti. Mentre le densità calcolate analiticamente possono richiedere trasformazioni matematiche complesse, i campioni MCMC possono essere facilmente trasformati con operazioni dirette.
In conclusione, il campionamento MCMC non è solo una necessità quando le soluzioni analitiche sono introvabili, ma offre vantaggi in termini di facilità di manipolazione, flessibilità computazionale e applicabilità pratica.
56.11 Caso Normale-Normale con Soluzione Analitica
Applichiamo ora l’algoritmo di Metropolis al caso Normale-Normale di cui conosciamo la soluzione analitica. In pratica, ci poniamo il problema di capire quale valore di
- un’idea iniziale (prior) che dice che
dovrebbe stare attorno a 30, con una certa incertezza (deviazione standard 5), - e abbiamo i dati osservati (
) che ci danno informazioni aggiuntive su dove si trova davvero .
Ma non conosciamo esattamente la distribuzione a posteriori di
Vogliamo costruire una “nuvola” di valori plausibili per
Step 1. Partiamo da un punto.
x_prev <- xinit
-
xinit è il valore iniziale: il nostro “primo sospetto” su dove si trovi
. - È come partire da un punto sulla mappa (“Penso che
sia circa qui”).
Step 2. Proponiamo un nuovo punto vicino.
x_star <- rnorm(1, mean = x_prev, sd = 0.5)
- Immaginiamo di essere bendati e di provare a fare un piccolo passo a caso partendo da dove siamo ora.
- Quel passo è generato con una distribuzione normale centrata su x_prev e con una deviazione standard piccola (0.5): piccoli passi casuali attorno al punto attuale.
Nota intuitiva:
Il valore 0.5 decide quanto “grandi” o “piccoli” sono i nostri passi. Più è grande, più possiamo saltare lontano; più è piccolo, più restiamo vicino.
Step 3. Calcoliamo quanto è “buono” il nuovo punto.
posterior(x_star, data)
posterior(x_prev, data)
- Ogni punto sulla mappa (
) ha un certo valore di plausibilità: quanto è probabile dati i dati osservati e il prior. - posterior(x_star, data) ci dice: “quanto è buono il nuovo punto?”
- posterior(x_prev, data) ci dice: “quanto era buono quello vecchio?”
Step 4. Decidiamo se accettare il nuovo punto.
Qui applichiamo il meccanismo di base dell’algoritmo di Metropolis per decidere sull’accettazione di un nuovo punto:
-
se il nuovo punto è migliore (cioè, la probabilità a posteriori è maggiore), allora lo accettiamo sicuramente (
, quindi ); -
se il nuovo punto è peggiore, possiamo comunque accettarlo con una certa probabilità:
- la probabilità di accettazione diminuisce all’aumentare di quanto il punto è “peggiore”;
- ciò è essenziale per non rimanere bloccati nei massimi locali.
In parole semplici:
- se troviamo un posto migliore, ci andiamo;
- se troviamo un posto peggiore, possiamo comunque andarci… ma tirando una monetina.
Step 5. Registriamo il punto attuale
Dopo aver deciso se accettare o meno il nuovo valore proposto, salviamo sempre un punto nella catena.
Ma attenzione:
- se la proposta è stata accettata, ci spostiamo al nuovo punto e lo registriamo;
- se la proposta è stata rifiutata, restiamo fermi e registriamo di nuovo la posizione attuale.
if (runif(1) < min(1, posterior(x_star, data) / posterior(x_prev, data))) {
x_prev <- x_star # accettiamo: ci spostiamo
}
samples[i] <- x_prev # salviamo dove ci troviamo ORA
In entrambi i casi, samples[i]
tiene traccia della posizione in cui ci troviamo dopo l’iterazione.
Intuizione: l’escursionista bendato.
Immaginiamo un’escursionista bendato che vuole esplorare un paesaggio fatto di colline di plausibilità (la distribuzione a posteriori):
- a ogni passo, prova a fare un salto in una nuova direzione (
x_star
); - se quel punto è più alto o non troppo peggiore, accetta di andarci e si sposta.
- se il punto è troppo brutto, rimane fermo dov’è;
- in ogni caso, segna nel diario la sua posizione attuale.
Ecco perché, quando guardiamo la catena, possiamo trovare valori ripetuti consecutivi: l’escursionista non si è mosso.
Questa caratteristica – il fatto che i campioni non siano tutti diversi – non è un errore, ma una proprietà fondamentale dell’algoritmo Metropolis: i campioni sono dipendenti e possono ripetersi.
Step 6. Ripetiamo tante volte.
for (i in seq_len(nsamp)) { ... }
- Più a lungo ripetiamo il processo (più iterazioni), più densa e accurata sarà la nostra approssimazione della distribuzione a posteriori di
. - Dopo un po’, i valori salvati formeranno un disegno della distribuzione plausibile di
.
** Riassunto in 3 frasi:**
- Partiamo da un valore sospettato di
.
- Partiamo da un valore sospettato di
- Facciamo piccoli passi casuali e decidiamo se accettarli in base a quanto sono “buoni” rispetto ai dati + prior.
- Dopo molti passi, la sequenza dei punti disegna la distribuzione a posteriori di
.
- Dopo molti passi, la sequenza dei punti disegna la distribuzione a posteriori di
** Dopo il sampling:**
- possiamo calcolare la media dei campioni = stima puntuale di
; - possiamo costruire un intervallo di credibilità = incertezza su
; - possiamo disegnare un istogramma dei campioni = forma della distribuzione a posteriori.
Anche se l’algoritmo di Metropolis può sembrare “rozzo” (tanti piccoli passi + accettare/rifiutare), funziona benissimo ed è uno dei motivi per cui oggi possiamo applicare la statistica bayesiana a modelli anche molto complessi.
Applichiamo dunque l’algoritmo di Metropolis all’esercizio in discussione. Iniziamo a definire le funzioni per il prior, la verosimiglianza e il posterior non normalizzato.
# Prior: Normal(30, 5^2)
prior <- function(mu) {
dnorm(mu, mean = 30, sd = 5)
}
# Likelihood: Normal(mu, sigma^2) con sigma calcolata dai dati
likelihood <- function(mu, data) {
sigma <- sd(data) # Deviazione standard dei dati
prod(dnorm(data, mean = mu, sd = sigma))
}
# Posterior non normalizzato
posterior <- function(mu, data) {
likelihood(mu, data) * prior(mu)
}
Implementiamo l’algoritmo di Metropolis per il caso normale-normale:
# Algoritmo di Metropolis
metropolis_for_normal <- function(nsamp, xinit, data) {
samples <- numeric(nsamp)
x_prev <- xinit
for (i in seq_len(nsamp)) {
x_star <- rnorm(1, mean = x_prev, sd = 0.5) # Proposta
if (runif(1) < min(1, posterior(x_star, data) / posterior(x_prev, data))) {
x_prev <- x_star
}
samples[i] <- x_prev
}
samples
}
Utilizziamo un campione di 30 valori BDI-II forniti da Zetsche et al. (2019):
# Dati osservati
y <- c(
26, 35, 30, 25, 44, 30, 33, 43, 22, 43, 24, 19, 39, 31, 25,
28, 35, 30, 26, 31, 41, 36, 26, 35, 33, 28, 27, 34, 27, 22
)
Esecuzione dell’algoritmo:
samples <- metropolis_for_normal(100000, mean(y), y)
Nel caso normale-normale, il posterior può essere calcolato analiticamente come segue:
# Parametri del prior
mu_prior <- 30
std_prior <- 5
var_prior <- std_prior^2
# Calcolo dei parametri posterior
n <- length(y)
sum_y <- sum(y)
var_data <- var(y)
mu_post <- (mu_prior / var_prior + sum_y / var_data) / (1 / var_prior + n / var_data)
var_post <- 1 / (1 / var_prior + n / var_data)
std_post <- sqrt(var_post)
mu_post
#> [1] 30.88
std_post
#> [1] 1.173
Visualizziamo i risultati con un istogramma dei campioni MCMC e la curva della distribuzione analitica:
# Assicuriamoci che palette_set1 sia disponibile
if (!exists("palette_set1")) {
palette_set1 <- RColorBrewer::brewer.pal(3, "Set1")
names(palette_set1) <- c("uno", "due", "tre")
}
# Campioni post burn-in
burnin <- floor(length(samples) * 0.5)
post_samples <- samples[-seq_len(burnin)]
# Dati per la curva analitica
x <- seq(mu_post - 4 * std_post, mu_post + 4 * std_post, length.out = 1000)
analytical_posterior <- dnorm(x, mean = mu_post, sd = std_post)
# Colori coerenti con la palette
colore_hist <- unname(palette_set1[[1]]) # es. blu Set1
colore_line <- unname(palette_set1[[2]]) # es. rosso Set1
# Creazione del grafico
ggplot() +
geom_histogram(
aes(x = post_samples, y = after_stat(density)),
bins = 30,
fill = colore_hist,
alpha = 0.4
) +
geom_line(
aes(x = x, y = analytical_posterior),
color = colore_line,
linewidth = 1
) +
labs(
title = "Distribuzione a posteriori: MCMC vs analitico",
x = expression(mu),
y = "Densità"
) +
theme(
plot.title = element_text(hjust = 0.5),
legend.position = "top"
)
Troviamo le proprietà del Posterior derivato con MCMC:
mean(samples)
#> [1] 30.91
sd(samples)
#> [1] 1.177
In conclusione, questo esempio illustra l’applicazione dell’algoritmo di Metropolis per la stima di una distribuzione a posteriori nel caso Normale-Normale e dimostra come confrontare i risultati del campionamento con la soluzione analitica, confermando così la coerenza tra le due approcci.
56.12 Riflessioni Conclusive
In molti casi, la distribuzione a posteriori dei parametri di un modello statistico non ha una forma analitica risolvibile. Per affrontare questa limitazione, si utilizzano metodi Monte Carlo basati su catene di Markov (MCMC). Questi algoritmi permettono di campionare efficacemente dalla distribuzione a posteriori, anche per modelli complessi, generando una sequenza di valori che approssima la distribuzione desiderata. L’algoritmo di Metropolis-Hastings (Hastings, 1970), un’estensione dell’algoritmo di Metropolis originale (Metropolis et al., 1953), è uno dei metodi MCMC più ampiamente utilizzati.
In sintesi, l’algoritmo segue questi passaggi principali:
- Generazione del nuovo stato proposto: Si crea un nuovo stato vicino a quello corrente utilizzando una distribuzione di proposta.
- Confronto tra densità posteriori: Si confrontano le densità a posteriori del nuovo stato proposto e dello stato corrente.
- Accettazione probabilistica: Il nuovo stato viene sempre accettato se ha una densità posteriore maggiore, oppure accettato con una certa probabilità se ha una densità minore.
- Burn-in e tasso di accettazione: I primi campioni vengono scartati (fase di burn-in) per garantire che la catena abbia raggiunto la distribuzione stazionaria, e si monitora il tasso di accettazione per ottimizzare l’efficienza del campionamento.
Questo approccio consente di ottenere campioni che approssimano la distribuzione a posteriori, ma l’algoritmo di Metropolis può presentare limiti di efficienza, soprattutto per problemi ad alta dimensionalità o distribuzioni con geometrie complesse. Un aspetto cruciale è il tasso di accettazione, che rappresenta il rapporto tra il numero di proposte accettate e il numero totale di proposte. Un tasso troppo basso può indicare che la catena esplora lo spazio dei parametri in modo inefficiente, mentre un tasso troppo alto può segnalare che i passi effettuati sono troppo piccoli per consentire una buona esplorazione.
Rispetto alle varianti più moderne, l’algoritmo di Metropolis tende a essere meno efficiente. Metodi come il No-U-Turn Sampler (NUTS) e l’Hamiltonian Monte Carlo (HMC) offrono importanti miglioramenti, specialmente in spazi di parametri di grandi dimensioni. NUTS, ad esempio, viene utilizzato in strumenti avanzati come Stan e PyMC (Hoffman et al., 2014), permettendo un’esplorazione più rapida e accurata della distribuzione a posteriori.
Tra gli altri algoritmi MCMC degni di nota troviamo il campionatore di Gibbs (Geman & Geman, 1984) e l’Hamiltonian Monte Carlo (Duane et al., 1987). Questi metodi, insieme a Metropolis-Hastings, formano la base di numerose tecniche moderne per il campionamento da distribuzioni complesse. Per un approfondimento dettagliato sulle tecniche MCMC, si consiglia di consultare Hanada & Matsuura (2022).
Informazioni sull’Ambiente di Sviluppo
sessionInfo()
#> R version 4.5.1 (2025-06-13)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS Sequoia 15.6
#>
#> Matrix products: default
#> BLAS: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.1
#>
#> locale:
#> [1] C/UTF-8/C/C/C/C
#>
#> time zone: Europe/Zagreb
#> tzcode source: internal
#>
#> attached base packages:
#> [1] stats graphics grDevices utils datasets methods base
#>
#> other attached packages:
#> [1] reshape2_1.4.4 cmdstanr_0.9.0 pillar_1.11.0
#> [4] tinytable_0.11.0 patchwork_1.3.1 ggdist_3.3.3
#> [7] tidybayes_3.0.7 bayesplot_1.13.0 ggplot2_3.5.2
#> [10] reliabilitydiag_0.2.1 priorsense_1.1.0 posterior_1.6.1
#> [13] loo_2.8.0 rstan_2.32.7 StanHeaders_2.32.10
#> [16] brms_2.22.0 Rcpp_1.1.0 conflicted_1.2.0
#> [19] janitor_2.2.1 matrixStats_1.5.0 modelr_0.1.11
#> [22] tibble_3.3.0 dplyr_1.1.4 tidyr_1.3.1
#> [25] rio_1.2.3 here_1.0.1
#>
#> loaded via a namespace (and not attached):
#> [1] gridExtra_2.3 inline_0.3.21 sandwich_3.1-1
#> [4] rlang_1.1.6 magrittr_2.0.3 multcomp_1.4-28
#> [7] snakecase_0.11.1 compiler_4.5.1 vctrs_0.6.5
#> [10] stringr_1.5.1 pkgconfig_2.0.3 arrayhelpers_1.1-0
#> [13] fastmap_1.2.0 backports_1.5.0 labeling_0.4.3
#> [16] rmarkdown_2.29 ps_1.9.1 purrr_1.1.0
#> [19] xfun_0.52 cachem_1.1.0 jsonlite_2.0.0
#> [22] broom_1.0.9 parallel_4.5.1 R6_2.6.1
#> [25] stringi_1.8.7 RColorBrewer_1.1-3 lubridate_1.9.4
#> [28] estimability_1.5.1 knitr_1.50 zoo_1.8-14
#> [31] pacman_0.5.1 R.utils_2.13.0 Matrix_1.7-3
#> [34] splines_4.5.1 timechange_0.3.0 tidyselect_1.2.1
#> [37] abind_1.4-8 yaml_2.3.10 codetools_0.2-20
#> [40] curl_6.4.0 processx_3.8.6 pkgbuild_1.4.8
#> [43] lattice_0.22-7 plyr_1.8.9 withr_3.0.2
#> [46] bridgesampling_1.1-2 coda_0.19-4.1 evaluate_1.0.4
#> [49] survival_3.8-3 RcppParallel_5.1.10 tensorA_0.36.2.1
#> [52] checkmate_2.3.2 stats4_4.5.1 distributional_0.5.0
#> [55] generics_0.1.4 rprojroot_2.1.0 rstantools_2.4.0
#> [58] scales_1.4.0 xtable_1.8-4 glue_1.8.0
#> [61] emmeans_1.11.2 tools_4.5.1 data.table_1.17.8
#> [64] mvtnorm_1.3-3 grid_4.5.1 QuickJSR_1.8.0
#> [67] colorspace_2.1-1 nlme_3.1-168 cli_3.6.5
#> [70] svUnit_1.0.6 Brobdingnag_1.2-9 V8_6.0.5
#> [73] gtable_0.3.6 R.methodsS3_1.8.2 digest_0.6.37
#> [76] TH.data_1.1-3 htmlwidgets_1.6.4 farver_2.1.2
#> [79] memoise_2.0.1 htmltools_0.5.8.1 R.oo_1.27.1
#> [82] lifecycle_1.0.4 MASS_7.3-65