12  Introduzione pratica a Stan

Introduzione

Nei capitoli precedenti abbiamo visto come l’algoritmo di Metropolis fornisca una soluzione generale al problema dell’inferenza bayesiana e come i linguaggi di programmazione probabilistica abbiano reso questa soluzione praticabile nella ricerca quotidiana. Ora è il momento di incontrare lo strumento che utilizzeremo concretamente nel nostro percorso: Stan.

Stan non è semplicemente un altro software statistico. È un linguaggio di programmazione probabilistica progettato per esprimere modelli bayesiani in modo chiaro e flessibile e per eseguire l’inferenza con algoritmi di campionamento allo stato dell’arte. In particolare, Stan utilizza varianti avanzate dell’Hamiltonian Monte Carlo (HMC), come l’algoritmo NUTS, che offrono efficienza e affidabilità molto superiori rispetto al semplice Metropolis.

Dal punto di vista concettuale, però, nulla cambia: la logica rimane la stessa che abbiamo già compreso. Stan non fa “magia”, ma implementa con grande efficacia ciò che Metropolis aveva già reso possibile. Per questo è importante vederlo come la naturale evoluzione pratica del percorso che abbiamo seguito fin qui.

Per la ricerca psicologica, Stan ha un vantaggio particolare. Molti dei modelli che ci interessano – modelli gerarchici, modelli di apprendimento, modelli dinamici – richiedono più parametri e strutture complesse. Implementarli a mano sarebbe quasi impossibile. Con Stan, invece, possiamo concentrarci sul modello teorico e tradurlo in codice relativamente semplice, lasciando al software la gestione dei dettagli computazionali.

Nei prossimi capitoli vedremo come iniziare a scrivere modelli in Stan, partendo da esempi elementari per arrivare a strutture più articolate. L’obiettivo non è soltanto imparare un nuovo linguaggio, ma acquisire la capacità di formalizzare i modelli psicologici come processi generativi, traducendoli in analisi statistiche riproducibili e trasparenti.

Panoramica del capitolo

• I blocchi del codice Stan.
  
• Scrivere e stimare modelli semplici con cmdstanr.
  
• Interpretare i risultati MCMC tramite riassunti e diagnostiche.
  
• Valutare la coerenza del modello con prior e posterior predictive check.

AttenzionePreparazione del Notebook

here::here("code", "_common.R") |> 
  source()
if (!requireNamespace("pacman")) install.packages("pacman")
pacman::p_load(cmdstanr, posterior, insight)

12.1 Programmazione probabilistica con Stan

Ogni programma Stan è organizzato in blocchi distinti, che corrispondono a funzioni precise (Nicenboim et al., 2025).

• Nel blocco data dichiariamo le variabili osservate che passiamo dall’esterno;
• in parameters indichiamo le quantità ignote che vogliamo stimare;
• nel blocco model specifichiamo le assunzioni probabilistiche — cioè le distribuzioni a priori e la verosimiglianza;
• in generated quantities possiamo calcolare misure derivate, come predizioni o log-likelihood, che non influiscono sulla stima ma sono utili per l’analisi successiva.

Possiamo pensare a un programma Stan come a una “ricetta”. Nel blocco data mettiamo gli ingredienti che già conosciamo (i dati osservati), in parameters dichiariamo gli ingredienti mancanti (i parametri da stimare), in model scriviamo le regole della preparazione (le distribuzioni a priori e la verosimiglianza) e in generated quantities prepariamo i contorni (diagnostiche, predizioni) che non cambiano la ricetta principale, ma la rendono più completa.

Questa architettura modulare rende Stan intuitivo e adattabile a un’ampia gamma di applicazioni statistiche.

12.1.1 Lavorare con Stan in R

L’interfaccia cmdstanr per R segue un workflow ben definito:

1. scrittura del modello in un file .stan,
   
2. compilazione del modello,
   
3. passaggio dei dati come lista R,
4. esecuzione del campionamento con sample(),
   
5. analisi dei risultati mediante pacchetti specializzati (posterior, bayesplot).

In pratica, lavorare con Stan da R segue sempre lo stesso schema. Prima si scrive il modello in un file .stan; poi lo si compila, cioè lo si traduce in un eseguibile; quindi si preparano i dati in una lista R con gli stessi nomi dichiarati nel modello; infine si lancia il campionamento con la funzione sample(). Una volta ottenuti i campioni a posteriori, possiamo analizzarli con pacchetti come posterior o bayesplot, che facilitano sia i riassunti numerici sia le visualizzazioni.

Stan utilizza un sistema di tipizzazione statica: tutte le variabili devono essere dichiarate con tipi specifici (int per interi, real per valori reali, vector per vettori) e possono includere vincoli (es. lower=0 per valori positivi). Questo approccio aumenta la robustezza del codice e previene errori comuni nella specificazione dei modelli.

12.2 Modello Beta–Binomiale

Come primo esempio costruiamo un modello molto semplice, ma già sufficiente per illustrare i principi fondamentali della programmazione in Stan. L’obiettivo è stimare la probabilità di successo \(\theta\) in una sequenza di prove Bernoulliane indipendenti.

Per rendere l’idea concreta, immaginiamo di lanciare un dado mille volte e di registrare il risultato come variabile dicotomica: assegniamo il valore 1 se esce il numero 6 (considerato “successo”), e 0 in tutti gli altri casi (“fallimento”). Se il dado fosse perfettamente equilibrato, la probabilità di successo sarebbe \(\theta = 1/6 \approx 0.167\). Tuttavia, nell’approccio bayesiano non assumiamo a priori che il dado sia equo: lasciamo che siano i dati, in combinazione con una distribuzione a priori esplicita, a informare il valore di \(\theta\).

Ecco un esempio di generazione dei dati in R:

set.seed(123)
n <- 1000
dice_df <- tibble(res = sample(1:6, size = n, replace = TRUE))
y <- as.integer(dice_df$res == 6)    # 1 se esce “6”, 0 altrimenti
mean(y)                              # frequenza relativa di “6”
#> [1] 0.164

Il modello statistico che proponiamo per rappresentare questa situazione è molto semplice:

• ogni osservazione segue una distribuzione Bernoulliana, \(y_i \sim \text{Bernoulli}(\theta)\) per \(i = 1, \dots, N\);
• ciò è equivalente a dire che il numero totale di successi \(k = \sum_i y_i\) segue una distribuzione Binomiale, \(k \sim \text{Binomiale}(N,\theta)\);
• come prior adottiamo una distribuzione uniforme su \([0,1]\), cioè \(\theta \sim \text{Beta}(1,1)\)¹.

12.2.1 Prima versione: modello Bernoulliano vettoriale

Il modo più diretto di tradurre il modello in Stan è scrivere la verosimiglianza come sequenza di esiti Bernoulliani. Questo approccio ha anche un valore didattico, perché mostra chiaramente la corrispondenza tra i dati osservati e il modello probabilistico. Stan, inoltre, vettorializza automaticamente le operazioni sugli array, rendendo il codice conciso:

stancode <- "
data {
  int<lower=1> N;                    // numero di prove
  array[N] int<lower=0, upper=1> y;  // esiti (0/1)
}
parameters {
  real<lower=0, upper=1> theta;      // probabilità di successo
}
model {
  theta ~ beta(1, 1);                // prior uniforme su [0,1]
  y ~ bernoulli(theta);              // verosimiglianza
}
"

Nel blocco data dichiariamo le variabili osservate (N e il vettore binario y); nel blocco parameters indichiamo il parametro da stimare, \(\theta\); infine, nel blocco model specifichiamo sia la distribuzione a priori sia la verosimiglianza.²

12.2.1.1 Perché array e non vector?

• In Stan, vector è usato solo per numeri reali (variabili continue a virgola mobile).
• Nel nostro caso gli esiti \(y\_i\) sono variabili discrete che assumono solo i valori 0 o 1. Per questo dobbiamo dichiararli come interi (int), e la struttura che raccoglie più interi in Stan è un array di interi.
• La sintassi array[N] int y; significa quindi “crea un array di lunghezza N, i cui elementi sono interi”.
• Non potremmo scrivere vector[N] y; perché un vector in Stan può contenere solo real, e non accetterebbe valori binari discreti come 0 o 1.

In altre parole, usiamo array ogni volta che i dati osservati sono categorici o interi (conteggi, esiti Bernoulli, ecc.), mentre vector è riservato a quantità continue come parametri, covariate normalizzate o variabili latenti.

12.2.2 Seconda versione: modello binomiale sui successi totali

L’approccio Bernoulliano visto in precedenza ha il pregio della chiarezza, ma può risultare ridondante: stiamo in realtà scrivendo la stessa formula mille volte, una per ciascun lancio del dado. In termini statistici, però, sappiamo che non è necessario conservare l’intera sequenza di zeri e uno: ai fini della stima di \(\theta\) conta solo il numero totale di successi osservati. Questa proprietà prende il nome di sufficienza della statistica \(k = \sum_i y_i\) per la distribuzione binomiale.

Se dunque nei mille lanci abbiamo osservato, ad esempio, 164 “6”, tutta l’informazione rilevante per stimare \(\theta\) è contenuta in quel singolo numero, non nella sequenza dettagliata dei lanci. La distribuzione binomiale ci permette di formalizzare questa idea in modo compatto, portando a una specificazione del modello più efficiente, ma del tutto equivalente sul piano inferenziale.

Ecco la traduzione in Stan:

stancode <- "
data {
  int<lower=1> N;                      // numero di prove
  array[N] int<lower=0, upper=1> y;    // sequenza di esiti (usata solo per calcolare k)
}
transformed data {
  int<lower=0, upper=N> k = sum(y);    // numero totale di successi
}
parameters {
  real<lower=0, upper=1> theta;        // probabilità di successo
}
model {
  theta ~ beta(1, 1);                  // prior uniforme
  k ~ binomial(N, theta);              // verosimiglianza sui successi totali
}
"

Questa versione concentra la verosimiglianza in un’unica riga, calcolata direttamente sul numero complessivo di successi. Le catene MCMC che otteniamo da questo modello coincidono, entro il rumore Monte Carlo, con quelle prodotte dalla versione Bernoulliana, ma richiedono meno operazioni e risultano quindi più efficienti dal punto di vista computazionale.

Come per la prima versione, possiamo arricchire il modello con un blocco generated quantities per ottenere log-likelihood o repliche predittive. Queste quantità derivate non modificano l’inferenza su \(\theta\), ma ci consentono di eseguire controlli diagnostici, confrontare modelli alternativi o visualizzare come il modello riproduce i dati osservati.

NotaPerché “Beta–Binomiale”?

Ricordiamo che, con prior \(\theta \sim \text{Beta}(a,b)\) e \(k \sim \text{Binomiale}(N,\theta)\), la posterior è

\[ \theta \mid y \sim \text{Beta}\big(a+k,\; b+N-k\big), \] con \(a=b=1\): \(\text{Beta}(1+k,\; 1+N-k)\).

Questo ci consente un controllo didattico: possiamo confrontare media e IC ottenuti da Stan con quelli della Beta a posteriori.

Esempio in R:

k <- sum(y)
a <- 1
b <- 1
post_mean_closed <- (a + k) / (a + b + n)
post_ci95_closed  <- qbeta(c(0.025, 0.975), a + k, b + n - k)
post_mean_closed; 
#> [1] 0.165
post_ci95_closed
#> [1] 0.142 0.188

Le stime via Stan (campioni MCMC della theta) devono coincidere (entro il rumore Monte Carlo) con queste quantità analitiche.

Una volta scritto il modello, il passo successivo è la compilazione. Stan traduce il codice in linguaggio C++ e lo trasforma in un piccolo eseguibile che potrà essere richiamato ogni volta che lanceremo le stime. Questo passaggio richiede qualche secondo solo la prima volta; in seguito, il modello compilato può essere riutilizzato con dataset diversi senza dover ricompilare da capo, con un notevole risparmio di tempo.

mod <- cmdstan_model(
  write_stan_file(stancode),
  compile = TRUE
)

Preparato l’eseguibile, dobbiamo passare a Stan i dati necessari. In questo caso servono due elementi: il numero totale di prove e il vettore con gli esiti dei lanci. È importante che i nomi e i tipi corrispondano esattamente a quanto dichiarato nel blocco data del modello Stan, altrimenti il programma restituirà un errore.

stan_data <- list(
  N = length(y),
  y = y
)
str(stan_data)
#> List of 2
#>  $ N: int 1000
#>  $ y: int [1:1000] 0 1 0 0 0 1 0 0 0 1 ...

A questo punto siamo pronti per il campionamento MCMC. Nella chiamata a sample() specifichiamo quante iterazioni dedicare alla fase di warmup (che serve per adattare l’algoritmo), quante iterazioni utilizzare effettivamente per l’inferenza, e quante catene indipendenti far partire in parallelo.³

fit1 <- mod$sample(
  data = stan_data,
  iter_warmup = 1000,
  iter_sampling = 2000,
  chains = 4,
  parallel_chains = 4,
  seed = 4790,
  refresh = 0
)
#> Running MCMC with 4 parallel chains...
#> 
#> Chain 1 finished in 0.0 seconds.
#> Chain 2 finished in 0.0 seconds.
#> Chain 3 finished in 0.0 seconds.
#> Chain 4 finished in 0.0 seconds.
#> 
#> All 4 chains finished successfully.
#> Mean chain execution time: 0.0 seconds.
#> Total execution time: 0.2 seconds.

Con il modello stimato, possiamo generare un riepilogo sintetico delle stime a posteriori. La funzione summary() di Stan mostra media, deviazione standard, quantili e diagnostiche come \(\hat R\) ed ESS.

fit1$summary(variables = "theta")
#> # A tibble: 1 × 10
#>   variable  mean median    sd   mad    q5   q95  rhat ess_bulk ess_tail
#>   <chr>    <dbl>  <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>    <dbl>    <dbl>
#> 1 theta    0.164  0.164 0.012 0.012 0.145 0.184 1.000 2443.381 3636.368

In alternativa, il pacchetto posterior permette di estrarre i campioni e calcolare con maggiore flessibilità le statistiche che ci interessano:

draws <- fit1$draws(variables = "theta", format = "draws_matrix")
posterior::summarise_draws(
  draws,
  mean, sd, ~quantile(.x, c(0.025, 0.5, 0.975)), rhat, ess_bulk, ess_tail
)
#> # A tibble: 1 × 9
#>   variable  mean    sd `2.5%` `50%` `97.5%`  rhat ess_bulk ess_tail
#>   <chr>    <dbl> <dbl>  <dbl> <dbl>   <dbl> <dbl>    <dbl>    <dbl>
#> 1 theta    0.164 0.012  0.142 0.164   0.188 1.000 2443.381 3636.368

Come regola generale, valori di \(\hat R\) molto vicini a 1 (idealmente < 1.01) e un numero elevato di campioni effettivi (ESS) indicano che le catene hanno esplorato bene la distribuzione a posteriori.

Oltre alle statistiche numeriche, le diagnostiche grafiche aiutano a valutare a colpo d’occhio la qualità del campionamento.

bayesplot::mcmc_trace(fit1$draws("theta"), n_warmup = 1000)

Nel traceplot, le catene dovrebbero mescolarsi bene senza mostrare trend persistenti; nel grafico delle densità, le distribuzioni delle diverse catene dovrebbero risultare sovrapposte.

In sintesi: la media a posteriori di \(\theta\) fornisce una stima puntuale della probabilità di ottenere un “6”, l’intervallo di credibilità quantifica l’incertezza residua, e le diagnostiche (R-hat vicino a 1, ESS alto, catene ben mescolate) garantiscono che i campioni siano affidabili per trarre conclusioni.

Confronto con la frequenza relativa

La frequenza mean(y) è solo una stima puntuale. L’analisi bayesiana restituisce un’intera distribuzione per \(\theta\), utile per propagare l’incertezza in previsioni e decisioni.

12.2.3 Distribuzione predittiva posteriore

Stimare \(\theta\) non è sufficiente: ciò che spesso ci interessa davvero è capire quali conseguenze pratiche derivano dai risultati. In altre parole, vogliamo sapere cosa il modello predice su dati futuri. Ad esempio: se rilanciassimo il dado altre 1000 volte, quanti “6” potremmo attenderci?

La risposta si ottiene tramite la distribuzione predittiva posteriore. Condizionatamente a un valore fissato di \(\theta\), il numero di successi nei nuovi lanci segue una distribuzione binomiale:

\[ y_{\text{rep}} \mid \theta \sim \text{Binomiale}(n_{\text{new}}, \theta). \]

Poiché non conosciamo \(\theta\) con certezza, ma solo la sua distribuzione a posteriori, la distribuzione predittiva si ottiene integrando la probabilità condizionata \(p(y\_{\text{rep}} \mid \theta)\) rispetto a \(p(\theta \mid y)\):

\[ p(y_{\text{rep}} \mid y) \;=\; \int p(y_{\text{rep}} \mid \theta)\, p(\theta \mid y)\, d\theta. \]

12.2.3.1 Procedura operativa

Il procedimento pratico è semplice e si articola in due passaggi:

1. Estrarre un valore da \(p(\theta \mid y)\). Supponiamo di ottenere \(\theta = 0.15\).

2. Simulare un dato replicato. Con quel valore generiamo un nuovo conteggio \(y\_{\text{rep}}\) da una binomiale, ad esempio con \(n\_{\text{new}}=1000\):

   rbinom(n = 1, size = 1000, prob = 0.15)
   #> [1] 160

Ripetendo il processo per molti valori di \(\theta\) estratti dalla distribuzione a posteriori, otteniamo un campione dalla distribuzione predittiva posteriore.

12.2.3.2 Codice R

# Numero di futuri lanci
n_new <- 1000

# Campioni da p(theta | y)
posterior_theta <- as.numeric(fit1$draws("theta", format = "matrix"))
posterior_theta[1:20]
#>  [1] 0.154 0.154 0.154 0.163 0.163 0.181 0.179 0.159 0.168 0.177 0.177 0.172
#> [13] 0.157 0.170 0.171 0.174 0.196 0.181 0.157 0.159

# Simulazione della distribuzione predittiva
yrep_count <- rbinom(
  n    = length(posterior_theta), 
  size = n_new, 
  prob = posterior_theta
)

# Dimensioni e prime osservazioni
length(yrep_count)
#> [1] 8000
head(yrep_count)
#> [1] 161 172 153 161 172 193

# Riassunti predittivi
mean_pred <- mean(yrep_count)  
ci_pred   <- quantile(yrep_count, c(0.025, 0.5, 0.975))

mean_pred
#> [1] 165
ci_pred
#>  2.5%   50% 97.5% 
#>   134   164   197

12.2.3.3 Visualizzazione

p025 <- ci_pred[1]; p50 <- ci_pred[2]; p975 <- ci_pred[3]

tibble(count = yrep_count) |>
  ggplot(aes(x = count)) +
  geom_histogram(bins = 30) +
  geom_vline(xintercept = mean_pred, linetype = "dashed") +
  geom_vline(xintercept = c(p025, p975), linetype = "dotted") +
  labs(
    x = "Numero di successi su n_new prove",
    y = "Frequenza",
    title = "Distribuzione predittiva posteriore",
    subtitle = "Linea tratteggiata = media predittiva; puntinate = intervallo al 95%"
  )

12.2.3.4 Interpretazione

L’istogramma mostra la variabilità attesa del numero di “6” su mille lanci futuri. La distribuzione è centrata intorno al valore atteso \(n\_{\text{new}} \cdot \mathbb{E}[\theta \mid y]\), cioè il numero medio di successi secondo la stima a posteriori.

Se il dado fosse perfettamente equo (\(\theta = 1/6\)), ci aspetteremmo circa 167 successi su 1000, con un intervallo predittivo attorno a 133–198. Un forte scostamento da questo scenario indicherebbe che il modello suggerisce un dado non equilibrato.

12.3 Dati continui: stima della media con varianza nota

Un problema molto comune nelle scienze psicologiche è la stima della media di popolazione di una variabile continua. Pensiamo, ad esempio, ai punteggi ottenuti a un test di intelligenza (QI). Per semplificare l’esposizione assumiamo che la deviazione standard della popolazione sia nota e pari a \(\sigma = 15\). Nella pratica questa informazione non è quasi mai disponibile, ma l’ipotesi ci permette di concentrare l’attenzione su un unico parametro incognito: la media \(\mu\). Ciò rende il modello più semplice e facilita il confronto tra l’impostazione classica e quella bayesiana.

Immaginiamo quindi di raccogliere i punteggi di \(n = 30\) persone, ipotizzando che la vera media sia 105:

set.seed(123)
n     <- 30
sigma <- 15
y_cont <- rnorm(n, mean = 105, sd = sigma)

tibble(y = y_cont) |>
  ggplot(aes(x = y)) +
  geom_histogram(bins = 15) +
  labs(x = "Punteggio", y = "Frequenza")

Il modello statistico che adottiamo ha due componenti fondamentali. La verosimiglianza descrive come si generano i dati: ogni osservazione \(y_i\) è distribuita secondo una Normale con media incognita \(\mu\) e deviazione standard nota \(\sigma\):

\[ y_i \sim \mathcal{N}(\mu, \sigma). \]

A questo si affianca il prior su \(\mu\), che formalizza le conoscenze disponibili prima di raccogliere i dati. Supponiamo che \(\mu\) sia distribuita normalmente attorno a \(\mu_0 = 100\), con deviazione standard \(\tau = 30\). Un prior così ampio riflette un’incertezza iniziale elevata, lasciando spazio ai dati per aggiornare in modo sostanziale le nostre convinzioni:

\[ \mu \sim \mathcal{N}(\mu_0, \tau). \]

Il modello combina dunque le informazioni fornite dai dati con quelle incorporate nel prior, producendo una distribuzione a posteriori per \(\mu\) che sintetizza entrambe le fonti di conoscenza.

La traduzione in Stan è la seguente:

stancode_norm <- "
data {
  int<lower=1> N;
  vector[N] y;             // dati osservati
  real<lower=0> sigma;     // deviazione standard nota
  real mu0;                // media del prior su mu
  real<lower=0> tau;       // sd del prior su mu
}
parameters {
  real mu;                 // media: parametro da stimare
}
model {
  mu ~ normal(mu0, tau);   // prior su mu
  y  ~ normal(mu, sigma);  // verosimiglianza
}
"

I dati devono essere preparati in R in un oggetto coerente con quanto richiesto dal blocco data:

stan_data2 <- list(
  N     = length(y_cont),
  y     = y_cont,
  sigma = sigma,
  mu0   = 100,
  tau   = 30
)

Dopo aver compilato il modello, si può avviare il campionamento MCMC:

mod2 <- cmdstanr::cmdstan_model(write_stan_file(stancode_norm), compile = TRUE)

fit2 <- mod2$sample(
  data = stan_data2,
  iter_warmup     = 1000,
  iter_sampling   = 4000,
  chains          = 4,
  parallel_chains = 4,
  seed            = 4790,
  refresh         = 0
)

I risultati possono essere riassunti in forma compatta:

print(fit2$summary(variables = "mu"), n = Inf)
#> # A tibble: 1 × 10
#>   variable    mean  median    sd   mad     q5     q95  rhat ess_bulk ess_tail
#>   <chr>      <dbl>   <dbl> <dbl> <dbl>  <dbl>   <dbl> <dbl>    <dbl>    <dbl>
#> 1 mu       104.240 104.236 2.721 2.740 99.751 108.710 1.001 5499.842 7540.041

oppure, con il pacchetto posterior:

draws2 <- fit2$draws(variables = "mu", format = "draws_matrix")
posterior::summarise_draws(
  draws2, mean, sd, ~quantile(.x, c(0.025, 0.5, 0.975)), rhat, ess_bulk, ess_tail
)
#> # A tibble: 1 × 9
#>   variable    mean    sd `2.5%`   `50%` `97.5%`  rhat ess_bulk ess_tail
#>   <chr>      <dbl> <dbl>  <dbl>   <dbl>   <dbl> <dbl>    <dbl>    <dbl>
#> 1 mu       104.240 2.721 98.916 104.236 109.627 1.001 5499.842 7540.041

L’interpretazione è diretta. Valori di \(\hat R\) vicini a 1 e un effective sample size elevato indicano che le catene sono ben miscelate e le stime affidabili. La media a posteriori di \(\mu\) rappresenta la miglior stima puntuale della media di popolazione. L’intervallo di credibilità quantifica l’incertezza residua: nel nostro caso esso risulta piuttosto ampio, e ciò riflette la limitata numerosità campionaria (\(n = 30\)). Con un campione più grande l’intervallo si restringerebbe, perché disporremmo di maggiore informazione empirica.

In conclusione, l’analisi bayesiana non fornisce soltanto una stima centrale, ma restituisce l’intera distribuzione di plausibilità dei valori possibili per la media della popolazione, offrendo così un quadro completo dell’incertezza associata ai dati raccolti.

12.3.1 Controllo analitico della coniugatezza

Il modello che abbiamo specificato ha una proprietà molto comoda: la coniugatezza. Quando usiamo un prior normale per la media \(\mu\) e assumiamo nota la varianza \(\sigma^2\), anche la distribuzione a posteriori di \(\mu\) rimane normale. Questo ci permette di calcolare media e varianza della posterior in forma chiusa, senza bisogno di simulazioni MCMC. In particolare, la varianza e la media della distribuzione a posteriori si ottengono come:

\[ \text{Var}(\mu \mid y) \;=\; \left(\frac{n}{\sigma^2} + \frac{1}{\tau^2}\right)^{-1}, \quad \mathbb{E}[\mu \mid y] \;=\; \text{Var}(\mu \mid y)\,\left(\frac{n\bar y}{\sigma^2} + \frac{\mu_0}{\tau^2}\right). \]

Con i nostri dati possiamo calcolare questi valori direttamente:

ybar    <- mean(y_cont)
tau2    <- 30^2
sigma2  <- sigma^2
post_var  <- 1 / (n / sigma2 + 1 / tau2)
post_sd   <- sqrt(post_var)
post_mean <- post_var * (n * ybar / sigma2 + 100 / tau2)

c(post_mean_closed = post_mean, post_sd_closed = post_sd)
#> post_mean_closed   post_sd_closed 
#>           104.26             2.73

Ora possiamo confrontare questi risultati analitici con quanto ottenuto tramite campionamento MCMC:

mu_draws <- as.numeric(draws2[, "mu"])
c(mean_mcmc = mean(mu_draws), sd_mcmc = sd(mu_draws))
#> mean_mcmc   sd_mcmc 
#>    104.24      2.72

Le due soluzioni dovrebbero concordare molto bene: la formula chiusa ci dà il risultato “esatto”, mentre l’MCMC lo approssima tramite simulazione. Le piccole discrepanze sono attribuibili al normale rumore Monte Carlo, che diminuisce all’aumentare del numero di campioni.

12.3.2 Visualizzazione

Una volta stimato il modello, è utile osservare graficamente i campioni MCMC per verificare sia la qualità del campionamento sia la forma della distribuzione a posteriori. Il pacchetto bayesplot fornisce strumenti immediati per questo scopo. Un primo passo è guardare l’istogramma dei campioni:

bayesplot::mcmc_hist(fit2$draws("mu"))

Questo grafico mostra la distribuzione stimata della media \(\mu\): non un singolo numero, ma un ventaglio di valori plausibili con le loro frequenze relative.

In sintesi, questo esempio rappresenta una sorta di “test bayesiano della media” con deviazione standard nota e un prior normale ampio su \(\mu\). A differenza dell’approccio frequentista, non otteniamo soltanto una stima puntuale o un intervallo, ma un’intera distribuzione a posteriori della media. Questo è un vantaggio importante: l’incertezza stimata può essere comunicata in modo trasparente e, soprattutto, può essere propagata nelle fasi successive dell’analisi, ad esempio in previsioni o decisioni basate sul modello.

12.3.3 Intervalli di credibilità

Nell’inferenza bayesiana non otteniamo un singolo numero come stima del parametro, ma un’intera distribuzione a posteriori che descrive l’incertezza residua. Per riassumere questa distribuzione in modo interpretabile si ricorre agli intervalli di credibilità, che indicano con quale probabilità il parametro cade in un certo intervallo dato il modello e i dati osservati. Dire che l’intervallo di credibilità al 94% per \(\mu\) va, ad esempio, da 101 a 110 significa che, condizionatamente alle assunzioni del modello e ai dati raccolti, c’è una probabilità del 94% che la media si trovi in quell’intervallo. Questa interpretazione è immediata e naturale, molto diversa da quella degli intervalli di confidenza frequentisti, che si riferiscono invece alla frequenza di copertura in ipotetiche ripetizioni del campione.

Gli intervalli di credibilità possono essere definiti in modi diversi. L’Equal-Tailed Interval (ETI) corrisponde all’intervallo che lascia la stessa probabilità nelle due code della distribuzione: ad esempio, in un intervallo al 94%, il 3% delle probabilità rimane sotto il limite inferiore e il 3% sopra quello superiore. Questo intervallo ha il pregio di essere stabile rispetto a trasformazioni monotone del parametro, ma può risultare inefficiente se la distribuzione è fortemente asimmetrica, perché finisce per includere regioni a bassa densità. L’Highest Density Interval (HDI), invece, è definito come l’intervallo più corto che contiene la probabilità prefissata. L’HDI si adatta meglio a distribuzioni asimmetriche, perché si concentra nelle zone di maggiore densità, ma ha lo svantaggio di non essere invariante a trasformazioni monotone e, in distribuzioni multimodali, può addirittura frammentarsi in più sotto-intervalli. Nei casi più semplici, quando la distribuzione a posteriori è simmetrica e unimodale, come per una Normale, ETI e HDI coincidono.

12.3.3.1 Esempio pratico con i campioni di \(\mu\)

Per chiarire l’uso di questi strumenti possiamo considerare i campioni MCMC della media \(\mu\) ottenuti dal modello normale con varianza nota. Calcolando sia l’ETI sia l’HDI al 94% si osserva che i due intervalli risultano praticamente identici, perché la distribuzione è simmetrica.

mu_draws <- as.numeric(fit2$draws("mu"))

Calcoliamo sia ETI che HDI al 94% con il pacchetto bayestestR:

eti94 <- bayestestR::eti(mu_draws, ci = 0.94)
hdi94 <- bayestestR::hdi(mu_draws, ci = 0.94)

eti94
#> 94% ETI: [99.12, 109.43]
hdi94
#> 94% HDI: [99.22, 109.48]

12.3.3.2 Visualizzazione

Con bayesplot::mcmc_areas() possiamo visualizzare gli intervalli centrali (ETI) e aggiungere i limiti HDI come linee verticali:

p <- bayesplot::mcmc_areas(fit2$draws("mu"), prob = 0.94) +
  xlab(expression(mu)) + ylab("Densità")

p + 
  geom_vline(xintercept = hdi94$CI_low,  linetype = "dashed") +
  geom_vline(xintercept = hdi94$CI_high, linetype = "dashed")

12.3.3.3 Posterior asimmetrica: esempio Beta

In distribuzioni asimmetriche, invece, le differenze emergono chiaramente. Un esempio utile è la distribuzione Beta(6,2), che mostra una coda a sinistra: in questo caso l’HDI risulta più corto e aderente alla parte densa della distribuzione, mentre l’ETI si estende ulteriormente verso la coda meno informativa.

set.seed(123)
theta_draws <- rbeta(5000, shape1 = 6, shape2 = 2)

eti_beta <- bayestestR::eti(theta_draws, ci = 0.94)
hdi_beta <- bayestestR::hdi(theta_draws, ci = 0.94)

eti_beta
#> 94% ETI: [0.43, 0.96]
hdi_beta
#> 94% HDI: [0.49, 0.98]

ggplot(data.frame(theta = theta_draws), aes(x = theta)) +
  geom_density() +
  geom_vline(xintercept = c(eti_beta$CI_low, eti_beta$CI_high), linewidth = 0.7) +
  geom_vline(xintercept = c(hdi_beta$CI_low, hdi_beta$CI_high), 
             linewidth = 1.2, linetype = "dashed") +
  labs(x = expression(theta), y = "Densità")

12.3.3.4 Quale livello utilizzare? 89%, 94% o 95%?

Rimane infine la questione di quale livello di probabilità utilizzare. La tradizione frequentista ha reso familiare il 95%, e molti analisti continuano a preferirlo per ragioni di continuità. Altri autori, come McElreath, suggeriscono il 94%, mentre Kruschke propone l’89% come scelta più stabile con campioni ridotti. In realtà non esiste una soglia “corretta” in senso assoluto: ciò che conta è la coerenza interna dell’analisi e la chiarezza nel riportare sia il livello scelto sia la definizione di intervallo adottata. In presenza di distribuzioni asimmetriche può essere particolarmente istruttivo mostrare sia l’ETI sia l’HDI, perché offrono prospettive complementari sulla distribuzione a posteriori del parametro.

12.4 Test di ipotesi in prospettiva bayesiana

L’approccio bayesiano consente di formulare domande in termini diretti e intuitivi, ad esempio:

“Qual è la probabilità che la media superi una determinata soglia?”

Questa impostazione si discosta da quella frequentista, che restituisce soltanto un esito dicotomico di rifiuto o non rifiuto dell’ipotesi nulla. In ambito bayesiano, al contrario, si ottiene una misura continua della plausibilità di ciascuna ipotesi, che rende l’incertezza più trasparente e consente di esprimere giudizi graduati.

12.4.1 Probabilità a posteriori sopra una soglia

Riprendiamo l’esempio dei punteggi QI con modello normale a varianza nota. Una prima domanda semplice è calcolare

\[ P(\mu > 110 \mid y), \]

cioè la probabilità che la media di popolazione superi il valore 110.

Questa quantità si ricava facilmente dai campioni MCMC della distribuzione a posteriori di \(\mu\):

# Estraiamo i campioni MCMC per mu
mu_draws <- as.numeric(fit2$draws("mu"))

# Calcoliamo la probabilità che mu superi la soglia 110
p_gt_110 <- mean(mu_draws > 110)
p_gt_110
#> [1] 0.0193

Il risultato rappresenta la probabilità, condizionatamente a modello e dati, che \(\mu\) sia superiore a 110. Non si tratta dunque di un’affermazione sì/no, ma di una misura graduata della plausibilità.

Il valore può essere verificato anche in forma analitica, sfruttando la coniugatezza del modello normale-normale:

# Calcolo analitico con distribuzione normale coniugata
p_gt_110_closed <- 1 - pnorm(110, mean = post_mean, sd = post_sd)
c(MCMC = p_gt_110, ClosedForm = p_gt_110_closed)
#>       MCMC ClosedForm 
#>     0.0193     0.0176

I due risultati dovrebbero coincidere, salvo piccole discrepanze dovute al rumore Monte Carlo.

12.4.2 Decisione pratica e ROPE

La valutazione delle ipotesi acquista senso compiuto solo quando si traduce in una decisione. Una regola semplice consiste nell’agire come se \(\mu > 110\) qualora la probabilità a posteriori superi una soglia prefissata \(p^\*\), scelta in base al bilancio tra costi e benefici. Ad esempio, con \(p^\* = 0.90\) si decide di agire solo se \(P(\mu > 110 \mid y) \geq 0.90\), altrimenti si rinvia o si raccolgono nuovi dati. La soglia \(p^\*\) non è arbitraria: un valore troppo basso aumenta il rischio di un’azione ingiustificata, mentre un valore troppo alto espone al rischio di trascurare un’azione necessaria. L’approccio bayesiano ha il merito di rendere esplicito questo nesso tra evidenza statistica e conseguenze pratiche, permettendo di calibrare la decisione sulla base delle priorità del contesto.

Spesso, tuttavia, l’interesse non è stabilire se la media sia esattamente pari a un certo valore, ma se possa essere considerata praticamente equivalente a esso entro un margine di tolleranza ritenuto accettabile. Per questo si introduce la ROPE (Region of Practical Equivalence), che definisce un intervallo attorno al valore soglia. Nel nostro esempio poniamo:

\[ \text{ROPE} = [108,\,112]. \] L’idea non è più chiedersi se \(\mu = 110\), ma se \(\mu\) appartenga a un intervallo che, per il problema in esame, può essere considerato equivalente a 110.

12.4.3 Calcolo e interpretazione

A partire dai campioni MCMC possiamo stimare la probabilità che \(\mu\) cada al di sotto, all’interno o al di sopra della ROPE:

rope <- c(108, 112)

P_below <- mean(mu_draws <  rope[1])                       # Probabilità sotto ROPE
P_in    <- mean(mu_draws >= rope[1] & mu_draws <= rope[2]) # Probabilità dentro ROPE
P_above <- mean(mu_draws >  rope[2])                       # Probabilità sopra ROPE

c(P_below = P_below, P_in = P_in, P_above = P_above)
#> P_below    P_in P_above 
#> 0.91463 0.08250 0.00287

Se la distribuzione a posteriori è ben approssimata da una Normale, gli stessi valori possono essere calcolati anche in forma chiusa:

P_below_cf <- pnorm(rope[1], mean = post_mean, sd = post_sd)
P_in_cf    <- pnorm(rope[2], mean = post_mean, sd = post_sd) -
              pnorm(rope[1], mean = post_mean, sd = post_sd)
P_above_cf <- 1 - pnorm(rope[2], mean = post_mean, sd = post_sd)

c(P_below_cf = P_below_cf, P_in_cf = P_in_cf, P_above_cf = P_above_cf)
#> P_below_cf    P_in_cf P_above_cf 
#>    0.91498    0.08275    0.00226

Se la probabilità che \(\mu\) cada all’interno della ROPE è elevata, si conclude che la media è praticamente equivalente al valore soglia. Se invece la probabilità si concentra al di sopra o al di sotto dell’intervallo, si parla rispettivamente di un’evidenza che \(\mu\) sia superiore o inferiore in senso rilevante. Nel caso in cui nessuna delle tre categorie prevalga nettamente, l’esito è quello di incertezza: in tale situazione è opportuno raccogliere ulteriori dati o riconsiderare l’ampiezza della ROPE adottata.

12.4.4 Esempio di report

Con ROPE = \([108,112]\), si ottiene: \(P(\mu < 108) = 0.92\), \(P(108 \leq \mu \leq 112) = 0.08\), \(P(\mu > 112) = 0.003\). Questi risultati indicano che la media \(\mu\) è con alta probabilità inferiore a 110, nel senso pratico definito dalla ROPE.

12.5 Diagnostiche di campionamento

Dopo aver ottenuto i campioni dalla distribuzione a posteriori, il passo successivo è verificare la loro qualità. Un modello ben specificato e un campionamento MCMC efficiente producono catene che esplorano in modo completo e bilanciato lo spazio dei parametri. Per questo, nei riepiloghi forniti da Stan e dai pacchetti associati compaiono alcune statistiche diagnostiche fondamentali, che meritano attenzione. La quantità lp__ compare in ogni output: non è un parametro del modello, ma la log-densità a posteriori, utile come strumento di controllo.

fit2$summary(variables = c("mu","lp__"))
#> # A tibble: 2 × 10
#>   variable    mean  median    sd   mad      q5     q95  rhat ess_bulk ess_tail
#>   <chr>      <dbl>   <dbl> <dbl> <dbl>   <dbl>   <dbl> <dbl>    <dbl>    <dbl>
#> 1 mu       104.240 104.236 2.721 2.740  99.751 108.710 1.001 5499.842 7540.041
#> 2 lp__     -14.463 -14.195 0.705 0.313 -15.896 -13.967 1.001 7356.664 9175.822

Un primo indicatore chiave è la statistica di convergenza \(\hat{R}\) (o R-hat). Essa confronta la variabilità intra-catena con quella inter-catena: quando le catene hanno raggiunto lo stesso equilibrio, oscillano nelle stesse regioni della distribuzione e \(\hat{R}\) tende a 1. Nel nostro esempio, i valori sono 1.001 per mu e 1.001 per lp__, praticamente indistinguibili da 1. In generale, valori inferiori a 1.01 sono considerati ottimali, mentre valori superiori a 1.05 segnalano possibili problemi di convergenza.

Un secondo gruppo di statistiche riguarda la dimensione del campione effettivo (ESS, effective sample size). Poiché i campioni MCMC sono correlati, l’ESS misura a quanti campioni indipendenti essi siano “equivalenti”. Nel nostro caso, per mu l’ESS è pari a 5500 (bulk) e 7540 (tail), mentre per lp__ è rispettivamente 7357 e 9176. Si tratta di valori elevati, che assicurano stime stabili anche per i quantili estremi. A titolo illustrativo, l’errore Monte Carlo sulla media di mu, calcolato come \(\text{sd}/\sqrt{\text{ESS}}\), è circa 0.04: un valore trascurabile rispetto alla deviazione standard della posterior, pari a 2.72.

Oltre agli indici numerici, è utile considerare la forma della distribuzione stimata. Nel riepilogo, la media (104.24) e la mediana (104.24) coincidono quasi perfettamente, segnalando una distribuzione sostanzialmente simmetrica. La deviazione standard (2.72) e la MAD (2.74) risultano molto vicine, confermando che la variabilità è ben catturata. I quantili dal 5° al 95° definiscono un intervallo credibile al 90% compreso tra circa e : ne consegue che la soglia 110 si colloca al di sopra di questo intervallo, con bassa probabilità di essere superata.

Per un riscontro visivo della dinamica delle catene, il traceplot deve mostrare linee ben mescolate, senza trend o stazionamenti anomali:

bayesplot::mcmc_trace(fit2$draws("mu"), n_warmup = 1000)

Un’ulteriore verifica consiste nel confronto delle densità tra catene:

bayesplot::mcmc_dens_overlay(fit2$draws("mu"))

La buona sovrapposizione delle curve indica che le catene hanno raggiunto la stessa distribuzione stazionaria.

Infine, lp__ — la log-densità a posteriori — non è un parametro di interesse sostantivo ma uno strumento diagnostico: anche per lp__ gli indici \(\hat{R}\) e ESS risultano eccellenti (1.001, 7357/9176), confermando un’esplorazione regolare dello spazio dei parametri.

In sintesi, tutti gli indici diagnostici — \(\hat{R}\) prossimo a 1, ESS elevati, coerenza tra media, mediana e misure di dispersione, oltre a traceplot e densità sovrapposte — convergono verso la stessa conclusione: il campionamento MCMC è stato stabile ed efficiente. Possiamo quindi considerare le stime ottenute come una rappresentazione affidabile della distribuzione a posteriori prevista dal nostro modello.

12.6 Prior e Posterior Predictive Check

Un passaggio cruciale nell’analisi bayesiana consiste nel verificare la capacità del modello di generare dati plausibili. Questo controllo può avvenire in due momenti distinti.

Prima dell’osservazione dei dati si esegue il prior predictive check, che serve a valutare se il prior scelto sia compatibile con valori ragionevoli. Se, ad esempio, i priors producono simulazioni di punteggi di QI del tutto inverosimili, è segno che occorre rivederne la specificazione.

Dopo l’aggiornamento con i dati osservati, si procede al posterior predictive check. In questa fase l’attenzione si sposta sulla capacità del modello stimato di riprodurre l’evidenza empirica: i dati simulati a posteriori dovrebbero assomigliare ai dati reali, sia nelle tendenze centrali sia nella variabilità. Se emergono scostamenti sistematici, ciò indica che il modello non riesce a catturare adeguatamente il fenomeno e richiede un perfezionamento.

Questi concetti sono già stati introdotti in precedenza; qui ci concentreremo sulla loro implementazione pratica utilizzando l’ecosistema di Stan.

12.6.1 Prior Predictive Check

Il primo passo, prima ancora di osservare i dati, è verificare se le nostre assunzioni a priori producono valori plausibili. In termini concreti: il prior scelto per \(\mu\) genera punteggi di QI che hanno senso rispetto al dominio del problema?

Il modello è definito come

\[ y_i \mid \mu \sim \mathcal{N}(\mu, \sigma), \qquad \mu \sim \mathcal{N}(\mu_0, \tau). \]

Combinando la verosimiglianza con il prior, otteniamo la distribuzione predittiva a priori per una singola osservazione:

\[ y_i \sim \mathcal{N}\!\Big(\mu_0,\; \sqrt{\sigma^2 + \tau^2}\Big). \] Questa distribuzione descrive quali valori ci aspetteremmo prima di osservare alcun dato. Se produce valori estremamente inverosimili (ad esempio QI < 40 o > 160), significa che il prior è troppo ampio o mal centrato. Se invece concentra quasi tutta la massa in un intervallo ristretto, rischia di essere eccessivamente informativo, lasciando poco spazio ai dati.

12.6.1.1 Implementazione in Stan

Per eseguire un prior predictive check possiamo utilizzare lo stesso modello Stan destinato all’inferenza, con una piccola modifica. Introduciamo una variabile booleana compute_likelihood che decide se attivare o meno la riga y ~ normal(mu, sigma);. In questo modo possiamo “spegnere” la verosimiglianza e simulare esclusivamente dal prior. Inoltre, nel blocco generated quantities generiamo delle repliche \(y_{\text{rep}}\) da confrontare con i dati osservati.

stancode_norm_ppc <- "
data {
  int<lower=0> N;
  vector[N] y;                 // usato solo se compute_likelihood=1
  real<lower=0> sigma;         // sd nota
  real mu0;                    // media del prior su mu
  real<lower=0> mu_prior_sd;   // sd del prior
  int<lower=0, upper=1> compute_likelihood; // 1 = attiva likelihood
}
parameters {
  real mu;
}
model {
  mu ~ normal(mu0, mu_prior_sd);
  if (compute_likelihood == 1) {
    y ~ normal(mu, sigma);
  }
}
generated quantities {
  vector[N] y_rep;
  vector[N] log_lik; 

  for (n in 1:N) {
    y_rep[n] = normal_rng(mu, sigma); 
    log_lik[n] = normal_lpdf(y[n] | mu, sigma); 
  }
}
"
stanmod_ppc <- cmdstan_model(write_stan_file(stancode_norm_ppc), compile = TRUE)

Per un controllo puro del prior impostiamo compute_likelihood = 0:

N_ppc <- length(y_cont)

stan_data_prior <- list(
  N = N_ppc,
  y = rep(0, N_ppc),   # placeholder, ignorato se compute_likelihood=0
  sigma = sigma,
  mu0 = 100,
  mu_prior_sd = 30,
  compute_likelihood = 0
)

fit_prior <- stanmod_ppc$sample(
  data = stan_data_prior,
  iter_warmup = 500,
  iter_sampling = 2000,
  chains = 4,
  parallel_chains = 4,
  seed = 4790,
  refresh = 500
)

Dopo il campionamento estraiamo le repliche \(y_{\text{rep}}\):

# Estrazione delle repliche
yrep_mat_prior <- posterior::as_draws_matrix(fit_prior$draws("y_rep"))

# Selezione delle colonne y_rep[1],...,y_rep[N]
yrep_mat_prior <- as.matrix(
  yrep_mat_prior[, paste0("y_rep[", 1:N_ppc, "]")]
)

Possiamo ora confrontare i dati osservati con alcune repliche simulate:

idx <- sample(seq_len(nrow(yrep_mat_prior)), 100)

bayesplot::ppc_dens_overlay(
  y = y_cont,
  yrep = yrep_mat_prior[idx, , drop = FALSE]
) 

12.6.1.2 Interpretazione

Se le distribuzioni simulate dal prior coprono la gamma dei dati reali, le nostre assunzioni iniziali sono plausibili. Se invece risultano troppo ampie o troppo strette, occorre ripensare la scelta di mu_prior_sd o la posizione del prior. In questo senso, il prior predictive check rappresenta un primo banco di prova fondamentale: ci permette di verificare se il modello, ancora prima di vedere i dati, è coerente con le conoscenze di base sul fenomeno studiato.

Nel grafico vediamo la distribuzione osservata dei punteggi (linea scura) confrontata con un insieme di repliche simulate dal prior (linee chiare). L’insieme delle simulazioni copre bene la variabilità dei dati reali: i punteggi osservati cadono in una regione ampiamente compatibile con quanto previsto a priori. Questo indica che il prior scelto è sufficientemente ampio da includere valori plausibili, senza però generare in prevalenza esiti palesemente irrealistici.

L’esercizio conferma quindi che le assunzioni iniziali sul parametro \(\mu\) non introducono vincoli eccessivi né distorsioni forti. Il modello, ancora prima di essere aggiornato con i dati, è coerente con la conoscenza di base sul fenomeno considerato (i punteggi di QI), e può essere usato come punto di partenza per l’inferenza bayesiana.

12.6.2 Posterior Predictive Check

Dopo aver verificato la plausibilità del prior, il passo successivo è valutare se il modello, una volta aggiornato con i dati, riesca a riprodurre in modo realistico l’evidenza empirica. Questo è lo scopo del posterior predictive check. Abbiamo già incontrato questo strumento nel caso del modello beta-binomiale; qui lo applichiamo al modello normale con prior normale, sfruttando le funzioni di Stan.

Per procedere, riattiviamo la verosimiglianza (compute_likelihood = 1) e stimiamo la distribuzione a posteriori di \(\mu\). Nel blocco generated quantities Stan genera repliche \(y_{\text{rep}}\), ossia nuovi dataset simulati sotto l’ipotesi che il modello e i parametri stimati siano corretti. Il confronto fra i dati osservati e queste repliche permette di valutare direttamente la capacità predittiva del modello.

Il principio è semplice: il prior predictive check serve a valutare le assunzioni prima di vedere i dati, mentre il posterior predictive check indaga se, dopo l’aggiornamento, il modello sia effettivamente in grado di riprodurre il fenomeno osservato.

stan_data_post <- list(
  N = length(y_cont),
  y = y_cont,
  sigma = sigma,
  mu0 = 100,
  mu_prior_sd = 30,
  compute_likelihood = 1
)

Il campionamento da Stan:

fit_post <- stanmod_ppc$sample(
  data = stan_data_post,
  iter_warmup = 1000,
  iter_sampling = 10000,
  chains = 4,
  parallel_chains = 4,
  seed = 4790,
  refresh = 1000
)

Una volta estratte le repliche, possiamo confrontarle con i dati reali:

y_rep <- fit_post$draws("y_rep", format = "matrix")
ppc_dens_overlay(y = stan_data_post$y, yrep = y_rep[1:100, ])

Se le distribuzioni simulate coprono bene quella dei dati osservati, il modello descrive adeguatamente il fenomeno. Al contrario, scostamenti sistematici — ad esempio repliche con media troppo bassa o varianza troppo elevata — segnalano che il modello non cattura pienamente la realtà e dovrebbe essere rivisto.

Il collegamento con la fase di verifica dei prior è immediato: un prior troppo vago può generare simulazioni poco credibili già in fase predittiva a priori, mentre un prior troppo rigido rischia di imporre vincoli eccessivi. Nel posterior predictive check, tuttavia, ciò che conta è il risultato finale: dopo l’aggiornamento con i dati, le simulazioni devono riflettere con realismo l’evidenza empirica.

Nel grafico osserviamo la distribuzione empirica dei punteggi (linea scura) affiancata a un insieme di repliche simulate dalla distribuzione a posteriori (linee chiare). A differenza del controllo sul prior, qui il modello è stato aggiornato con i dati, e le simulazioni riflettono sia le assunzioni iniziali sia l’informazione empirica raccolta.

Il confronto mostra una buona corrispondenza: le repliche coprono con realismo la forma della distribuzione osservata, senza evidenziare scostamenti sistematici. In particolare, la media delle repliche si allinea bene a quella dei dati reali, e la variabilità simulata risulta compatibile con quella empirica. Questo conferma quanto già emerso dall’analisi numerica: la stima a posteriori della media \(\mu\) (circa 104) è coerente con il campione osservato e attribuisce bassa probabilità all’ipotesi che \(\mu\) superi la soglia di 110.

Nel complesso, il posterior predictive check suggerisce che il modello normale con prior scelto descrive adeguatamente il fenomeno considerato (i punteggi QI in questo esempio). Si tratta quindi di una base affidabile per trarre conclusioni inferenziali e, se necessario, collegare i risultati a decisioni pratiche.

Riflessioni conclusive

Stan rappresenta il punto di approdo naturale del percorso seguito fin qui. Dopo aver introdotto la logica dell’inferenza bayesiana con algoritmi semplici come Metropolis e aver visto come i linguaggi probabilistici abbiano reso più accessibile questa logica, Stan ci mette ora a disposizione uno strumento maturo, capace di tradurre queste idee in pratica su modelli reali e complessi.

Il suo valore non risiede soltanto nella potenza computazionale, ma soprattutto nella possibilità di spostare l’attenzione dal calcolo all’idea scientifica. Con Stan possiamo formulare ipotesi psicologiche in termini formali, esplicitarne le assunzioni e sottoporle a verifica empirica in modo trasparente e riproducibile. Il ricercatore non è più costretto a concentrarsi sui dettagli algoritmici, ma può dedicarsi a ciò che conta davvero: la sostanza teorica del modello e la sua capacità di descrivere i dati.

In questo senso, Stan non è semplicemente un software, ma un vero e proprio ambiente di ricerca che promuove chiarezza, rigore e cumulatività. Nei capitoli che seguono impareremo a utilizzarlo prima con esempi semplici, per poi estendere gradualmente l’analisi a modelli più articolati. L’obiettivo non è solo acquisire familiarità con la sintassi, ma soprattutto interiorizzare un modo di pensare: quello che fa della modellazione bayesiana uno strumento essenziale della psicologia scientifica contemporanea.

NotaInformazioni sull’ambiente di sviluppo

sessionInfo()
#> R version 4.5.1 (2025-06-13)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS Tahoe 26.0
#> 
#> 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/UTF-8/C.UTF-8/C/C.UTF-8/C.UTF-8
#> 
#> time zone: Europe/Rome
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] insight_1.4.2         cmdstanr_0.9.0        pillar_1.11.1        
#>  [4] tinytable_0.13.0      patchwork_1.3.2       ggdist_3.3.3         
#>  [7] tidybayes_3.0.7       bayesplot_1.14.0      ggplot2_4.0.0        
#> [10] reliabilitydiag_0.2.1 priorsense_1.1.1      posterior_1.6.1      
#> [13] loo_2.8.0             rstan_2.32.7          StanHeaders_2.32.10  
#> [16] brms_2.23.0           Rcpp_1.1.0            sessioninfo_1.2.3    
#> [19] conflicted_1.2.0      janitor_2.2.1         matrixStats_1.5.0    
#> [22] modelr_0.1.11         tibble_3.3.0          dplyr_1.1.4          
#> [25] tidyr_1.3.1           rio_1.2.3             here_1.0.2           
#> 
#> 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.4        multcomp_1.4-28      
#>  [7] snakecase_0.11.1      ggridges_0.5.7        compiler_4.5.1       
#> [10] reshape2_1.4.4        systemfonts_1.2.3     vctrs_0.6.5          
#> [13] stringr_1.5.2         pkgconfig_2.0.3       arrayhelpers_1.1-0   
#> [16] fastmap_1.2.0         backports_1.5.0       labeling_0.4.3       
#> [19] utf8_1.2.6            rmarkdown_2.29        ps_1.9.1             
#> [22] ragg_1.5.0            purrr_1.1.0           xfun_0.53            
#> [25] cachem_1.1.0          jsonlite_2.0.0        broom_1.0.10         
#> [28] parallel_4.5.1        R6_2.6.1              stringi_1.8.7        
#> [31] RColorBrewer_1.1-3    lubridate_1.9.4       estimability_1.5.1   
#> [34] knitr_1.50            zoo_1.8-14            pacman_0.5.1         
#> [37] Matrix_1.7-4          splines_4.5.1         timechange_0.3.0     
#> [40] tidyselect_1.2.1      abind_1.4-8           yaml_2.3.10          
#> [43] codetools_0.2-20      curl_7.0.0            processx_3.8.6       
#> [46] pkgbuild_1.4.8        plyr_1.8.9            lattice_0.22-7       
#> [49] bayestestR_0.17.0     withr_3.0.2           bridgesampling_1.1-2 
#> [52] S7_0.2.0              coda_0.19-4.1         evaluate_1.0.5       
#> [55] survival_3.8-3        RcppParallel_5.1.11-1 tensorA_0.36.2.1     
#> [58] checkmate_2.3.3       stats4_4.5.1          distributional_0.5.0 
#> [61] generics_0.1.4        rprojroot_2.1.1       rstantools_2.5.0     
#> [64] scales_1.4.0          xtable_1.8-4          glue_1.8.0           
#> [67] emmeans_1.11.2-8      tools_4.5.1           data.table_1.17.8    
#> [70] mvtnorm_1.3-3         grid_4.5.1            QuickJSR_1.8.0       
#> [73] colorspace_2.1-1      nlme_3.1-168          cli_3.6.5            
#> [76] textshaping_1.0.3     svUnit_1.0.8          Brobdingnag_1.2-9    
#> [79] V8_7.0.0              gtable_0.3.6          digest_0.6.37        
#> [82] TH.data_1.1-4         htmlwidgets_1.6.4     farver_2.1.2         
#> [85] memoise_2.0.1         htmltools_0.5.8.1     lifecycle_1.0.4      
#> [88] MASS_7.3-65

Bibliografia

Nicenboim, B., Schad, D. J., & Vasishth, S. (2025). Introduction to Bayesian data analysis for cognitive science. CRC Press.

---

1. Vale la pena notare che, se in Stan dichiariamo un parametro vincolato all’intervallo \([0,1]\) ma non specifichiamo alcun prior, il linguaggio assegna automaticamente proprio questa distribuzione uniforme, che corrisponde a una Beta(1,1).

2. Dal punto di vista interno, ogni riga del tipo x ~ distribuzione(...) aggiunge la log-densità (o log-massa) alla quantità interna target, che rappresenta la log-posterior. Se vogliamo essere più espliciti, possiamo scrivere il codice in forma equivalente:

   target += beta_lpdf(theta | 1, 1);
   target += bernoulli_lpmf(y | theta);

   Questa seconda scrittura, sebbene meno compatta, è particolarmente utile quando vogliamo costruire verosimiglianze personalizzate o introdurre modifiche non standard.

3. La funzione sample() ha anche altri parametri importanti per ottimizzare il campionamento. Tra questi, adapt_delta aumenta la precisione dell’algoritmo di campionamento, riducendo gli errori di divergenza, mentre max_treedepth controlla la profondità massima dell’albero di campionamento. Utilizzarli è fondamentale per ottenere risultati stabili e affidabili, specialmente quando si lavora con modelli complessi o con dati problematici.