44  Aggiornare le credenze su un parametro: dal prior alla posterior

“The probability of any event is the ratio between the value at which an expectation depending on the happening of the event ought to be computed, and the value of the thing expected upon its happening.”

— Thomas Bayes, An Essay towards Solving a Problem in the Doctrine of Chances

Introduzione

In questo capitolo esploreremo l’inferenza bayesiana attraverso il modello binomiale, particolarmente adatto per analizzare dati dicotomici (successo/fallimento, sì/no, corretto/errato). Inizieremo con un caso molto semplice: una distribuzione a priori discreta, dove il parametro \(\theta\) (la probabilità di successo) può assumere solo un numero limitato di valori predefiniti. Questo ci permetterà di vedere, passo dopo passo, come le nostre credenze su \(\theta\) cambiano quando arrivano nuovi dati. Successivamente passeremo al caso più generale delle distribuzioni a priori continue, molto più adatte a modellizzare problemi reali. In particolare useremo la distribuzione Beta, che si combina perfettamente con la distribuzione Binomiale.

L’obiettivo è fornire una comprensione intuitiva ma rigorosa del processo di aggiornamento bayesiano, mettendo in luce sia la logica di fondo sia le potenzialità applicative.

Panoramica del capitolo

  • applicare l’aggiornamento bayesiano per affinare credenze;
  • rappresentare distribuzioni a priori (discrete e continue);
  • calcolare la verosimiglianza e aggiornare la distribuzione a priori;
  • derivare e interpretare la distribuzione a posteriori;
  • usare il metodo a griglia per approssimare la distribuzione a posteriori;
  • applicare il modello binomiale per stimare probabilità e incertezze;
  • calcolare medie, mode e intervalli di credibilità;
  • utilizzare la distribuzione Beta come prior continuo.
  • Leggere il settimo capitolo del testo di Albert & Hu (2019).
here::here("code", "_common.R") |> 
  source()

# Load packages
if (!requireNamespace("pacman")) install.packages("pacman")
pacman::p_load(HDInterval)

44.1 Verosimiglianza binomiale

La distribuzione binomiale descrive il numero di successi \(y\) in \(n\) prove indipendenti, ciascuna con probabilità di successo \(\theta\):

\[ p(y \mid \theta) = \text{Bin}(y \mid n, \theta) = \binom{n}{y} \theta^y (1 - \theta)^{n-y}, \]

dove \(\theta\) rappresenta la probabilità di successo per singola prova, \(y\) è il numero osservato di successi e \(n\) è il numero totale di prove (fissato a priori).

44.2 Esempio psicologico: riconoscimento di parole emotive

Un ricercatore in psicologia cognitiva conduce uno studio per stimare la proporzione di adulti in grado di riconoscere correttamente parole con contenuto emotivo positivo (es. gioia, orgoglio, speranza) all’interno di un compito di memoria a lungo termine. Ogni partecipante riceve lo stesso test standardizzato e viene classificato come:

  • 1 = corretto riconoscimento (successo), oppure
  • 0 = mancato riconoscimento (fallimento).

In un campione di 30 soggetti indipendenti, 22 hanno riconosciuto correttamente le parole emozionanti.

In questo scenario:

  • ciascun soggetto è un’unità di osservazione indipendente;
  • ciascun soggetto può essere considerato come una prova di una variabile di Bernoulli con probabilità di successo \(\theta\), dove \(\theta\) rappresenta la proporzione vera nella popolazione di adulti in grado di riconoscere correttamente parole emotive.

Pertanto, il numero totale di successi in \(n = 30\) soggetti può essere modellato come:

\[ Y \sim \text{Binomiale}(n = 30, \theta), \]

dove \(Y = 22\) è il numero di successi osservati.

44.2.1 Obiettivo inferenziale

Il nostro scopo è stimare \(\theta\): la probabilità che un adulto, scelto a caso dalla popolazione, riconosca correttamente le parole emotive. Nel quadro bayesiano combiniamo:

  • la verosimiglianza (dati osservati: 22 su 30),
  • una distribuzione a priori (credenze iniziali su \(\theta\)),

ottenendo la distribuzione a posteriori, che rappresenta la nostra conoscenza aggiornata.

44.3 Metodo basato su griglia

Il metodo su griglia è un approccio semplice per calcolare la distribuzione a posteriori. È didatticamente utile perché rende espliciti tutti i passaggi.

Passaggi:

  • definire una griglia di valori plausibili per \(\theta\) (da 0 a 1),
  • calcolare la verosimiglianza per ogni valore di \(\theta\),
  • moltiplicare verosimiglianza e priori,
  • normalizzare i valori ottenuti in modo che la somma sia 1.

Il risultato è una distribuzione a posteriori che mostra come l’evidenza aggiorna le credenze iniziali.

44.4 Aggiornamento bayesiano con una distribuzione a priori discreta

44.4.1 Costruzione della distribuzione a priori

In assenza di informazioni specifiche, possiamo assumere che tutti i valori di \(\theta\) siano ugualmente plausibili. Per implementare concretamente questo approccio:

  • definiamo un insieme discreto di valori possibili per \(\theta\): {0, 0.1, 0.2, …, 1},
  • assegniamo a ciascun valore la stessa probabilità a priori: \(p(\theta) = 1/11 \approx 0.09\).

Questa scelta rappresenta uno stato di massima incertezza iniziale, dove nessun valore di \(\theta\) risulta a priori più plausibile di altri.

44.4.2 Aggiornamento con i dati

Abbiamo osservato \(y = 22\) successi su \(n = 30\). Per ogni valore \(\theta\) nella griglia:

  1. calcoliamo la verosimiglianza binomiale: \(p(y \mid \theta) = \theta^{22}(1-\theta)^8\),
  2. moltiplichiamo per la probabilità a priori,
  3. normalizziamo dividendo per la somma totale di tutti i prodotti ottenuti.

Il risultato è una distribuzione a posteriori discreta che mostra come l’osservazione di 22 successi su 30 prove modifica le nostre credenze iniziali su \(\theta\). I valori più vicini a 22/30 \(\approx\) 0.7 riceveranno una maggiore probabilità a posteriori.

44.4.3 Interpretazione

  • Prima dei dati, ogni valore di \(\theta\) era ugualmente plausibile.
  • Dopo i dati, valori come \(\theta = 0.7\) o \(0.75\) hanno alta probabilità a posteriori.
  • Valori estremi (\(0.2\), \(0.9\)) diventano poco plausibili.

44.4.4 Implementazione in R

Definizione della griglia:

theta <- seq(0, 1, by = 0.1)  # Griglia di valori da 0 a 1 con passo 0.1

Quando non abbiamo informazioni preliminari, usiamo una distribuzione uniforme:

priori_unif <- rep(1 / length(theta), length(theta))  # Probabilità uniformi

Visualizziamo questa distribuzione:

ggplot(data.frame(theta, prob = priori_unif), aes(x = theta, y = prob)) +
  geom_col(width = 0.08) +
  labs(title = "Distribuzione a Priori Uniforme",
       x = expression(theta),
       y = "Densità di probabilità")

Se invece riteniamo più probabili valori centrali di \(\theta\):

priori_inf <- c(
  0, 0.05, 0.05, 0.05, 0.175, 0.175, 0.175, 0.175, 0.05, 0.05, 0.05
)

Visualizzazione:

ggplot(data.frame(theta, prob = priori_inf), aes(x = theta, y = prob)) +
  geom_col(width = 0.08) +
  labs(title = "Distribuzione a Priori Informativa",
       x = expression(theta),
       y = "Densità di probabilità")

Verosimiglianza:

verosimiglianza <- dbinom(22, size = 30, prob = theta)

Visualizzazione:

ggplot(data.frame(theta, prob = verosimiglianza), aes(x = theta, y = prob)) +
  geom_col(width = 0.08) +
  labs(title = "Funzione di Verosimiglianza",
       x = expression(theta),
       y = "L(θ|dati)")

Calcolo della distribuzione a posteriori:

posteriori_non_norm <- priori_inf * verosimiglianza
posteriori <- posteriori_non_norm / sum(posteriori_non_norm)  # Normalizzazione

Visualizzazione:

ggplot(data.frame(theta, prob = posteriori), aes(x = theta, y = prob)) +
  geom_col(width = 0.08) +
  labs(title = "Distribuzione a Posteriori",
       x = expression(theta),
       y = "P(θ|dati)")

Statistiche descrittive:

media_post <- sum(theta * posteriori)
var_post <- sum(theta^2 * posteriori) - media_post^2
moda_post <- theta[which.max(posteriori)]

cat("Media a posteriori:", round(media_post, 3),
    "\nVarianza a posteriori:", round(var_post, 3),
    "\nModa a posteriori:", moda_post)
#> Media a posteriori: 0.689 
#> Varianza a posteriori: 0.005 
#> Moda a posteriori: 0.7

L’implementazione illustra tre caratteristiche essenziali dell’inferenza bayesiana. La funzione di verosimiglianza attribuisce maggiore densità di probabilità a valori di \(\theta\) compresi tra 0.6 e 0.8, in accordo con l’evidenza empirica dei 22 successi osservati su 30 prove. La distribuzione a priori contribuisce in modo determinante alla configurazione della distribuzione a posteriori risultante. Il processo di aggiornamento bayesiano integra in modo rigoroso queste due fonti informative mediante l’applicazione del teorema di Bayes.

44.5 Aggiornamento bayesiano con una distribuzione a priori continua

Un’estensione naturale è usare una distribuzione continua come priori. La più adatta nel caso di proporzioni è la Beta:

  • supporto: \([0,1]\) (come \(\theta\));
  • coniugata della Binomiale (la posteriori è ancora una Beta);
  • parametri: \(\text{Beta}(\alpha, \beta)\).

Esempi:

  • \(\text{Beta}(2,2)\): priori simmetrica e non troppo informativa.
  • \(\text{Beta}(2,5)\): priori che privilegia valori bassi di \(\theta\).

44.5.1 Implementazione in R

Calcoliamo la densità della distribuzione \(\text{Beta}(2, 2)\) su una griglia fine di valori di \(\theta\):

theta <- seq(0, 1, length.out = 1000)
prior_beta_2_2 <- dbeta(theta, 2, 2) / sum(dbeta(theta, 2, 2))

Visualizzazione:

ggplot(data.frame(theta, prior = prior_beta_2_2), aes(x = theta, y = prior)) +
  geom_line(linewidth = 1.2) +
  labs(title = "Distribuzione a Priori Beta(2, 2)", x = expression(theta), y = "Densità")

Continuiamo con l’esempio precedente, con 22 successi in 30 prove. La verosimiglianza associata a ciascun valore di \(\theta\) è calcolata come:

likelihood <- dbinom(22, size = 30, prob = theta)
likelihood <- likelihood / sum(likelihood)  # normalizzazione opzionale

Poiché il prior è continuo, otteniamo la distribuzione a posteriori moltiplicando punto a punto la densità a priori per la verosimiglianza, e normalizzando:

posterior_unnorm <- prior_beta_2_2 * likelihood
posterior <- posterior_unnorm / sum(posterior_unnorm)

Visualizziamo le tre curve:

df <- data.frame(theta, prior = prior_beta_2_2, likelihood, posterior)

df_long <- df |>
  pivot_longer(cols = c("prior", "likelihood", "posterior"),
               names_to = "Distribuzione", values_to = "Densità")

ggplot(df_long, aes(x = theta, y = Densità, color = Distribuzione)) +
  geom_line(size = 1.2) +
  labs(title = "Aggiornamento Bayesiano con Prior Continua",
       x = expression(theta), y = "Densità") +
  theme(legend.position = "bottom")

Consideriamo ora una distribuzione a priori non simmetrica, Beta(2, 5), per rappresentare credenze che privilegiano valori bassi di \(\theta\).

prior_2_5 <- dbeta(theta, 2, 5) / sum(dbeta(theta, 2, 5))

La distribuzione a posteriori si ottiene moltiplicando la distribuzione a priori per la verosimiglianza e normalizzando il risultato:

posterior <- (prior_2_5 * likelihood) / sum(prior_2_5 * likelihood)
# Uniamo tutti i dati in un dataframe per ggplot2
dat <- tibble(
  theta,
  prior_2_5,
  likelihood,
  posterior
)

# Preparazione dei dati per il plot
long_data <- dat |>
  pivot_longer(
    cols = c(prior_2_5, likelihood, posterior),
    names_to = "distribution",
    values_to = "density"
  )

# Grafico
long_data |>
  ggplot(aes(x = theta, y = density, color = distribution)) +
  geom_line(size = 1.2) +
  labs(
    title = "Distribuzioni Bayesiane",
    x = expression(theta),
    y = "Densità"
  ) +
  theme(legend.position = "bottom")

Calcoliamo alcune quantità descrittive utili:

posterior_mean <- sum(theta * posterior)
posterior_sd <- sqrt(sum((theta^2) * posterior) - posterior_mean^2)
posterior_mode <- theta[which.max(posterior)]
posterior_mean; posterior_sd; posterior_mode
#> [1] 0.6486
#> [1] 0.07744
#> [1] 0.6567

Per stimare intervalli di credibilità, possiamo campionare dalla distribuzione a posteriori:

samples <- sample(theta, size = 10000, replace = TRUE, prob = posterior)
quantile(samples, probs = c(0.03, 0.97))  # intervallo al 94%
#>     3%    97% 
#> 0.4955 0.7878

Se desideriamo calcolare l’intervallo di densità più alta (HPDI), possiamo utilizzare pacchetti aggiuntivi come HDInterval.

# Calcolo HPDI (richiede il pacchetto HDInterval)
hdi(samples, credMass = 0.94)
#>  lower  upper 
#> 0.5025 0.7928 
#> attr(,"credMass")
#> [1] 0.94

Interpretazione. L’aggiornamento bayesiano con una distribuzione a priori continua fornisce una stima aggiornata di \(\theta\) che tiene conto sia della distribuzione a priori (conoscenza pregressa) sia della verosimiglianza (dati osservati). Nel nostro esempio, la curva a posteriori risulta spostata verso destra rispetto al prior simmetrico \(\text{Beta}(2,2)\), riflettendo l’evidenza di 22 successi su 30 prove. In alternativa, utilizzando un prior asimmetrico come \(\text{Beta}(2,5)\), la distribuzione a posteriori mostra un compromesso tra la tendenza iniziale a credere in basse probabilità di successo e l’evidenza empirica più ottimista.

Il metodo basato su griglia si estende anche ad altri modelli, come il caso della distribuzione normale. Per esempio, possiamo usare il metodo su griglia per stimare la media \(\mu\) di una popolazione, assumendo che i dati seguano una distribuzione normale con deviazione standard nota. Anche in questo contesto, la procedura prevede di:

  1. definire una griglia di valori plausibili per il parametro \(\mu\),
  2. calcolare la verosimiglianza dei dati per ciascun valore della griglia,
  3. specificare una distribuzione a priori per \(\mu\),
  4. combinare prior e verosimiglianza punto per punto,
  5. normalizzare i risultati per ottenere una distribuzione a posteriori.

Esempio: Stima bayesiana della media del QI.

Immaginiamo di condurre uno studio su bambini ad alto potenziale cognitivo. Supponiamo che i loro punteggi di QI siano distribuiti secondo una normale con deviazione standard nota \(\sigma = 5\). La media \(\mu\) è incognita e rappresenta il parametro di interesse.

Dati osservati. Abbiamo raccolto un piccolo campione (10 bambini):

set.seed(123)
campione <- round(rnorm(10, mean = 130, sd = 5))
campione
#>  [1] 127 129 138 130 131 139 132 124 127 128

La verosimiglianza congiunta. La verosimiglianza esprime quanto un certo valore del parametro \(\mu\) renda plausibili i dati osservati.

Se osserviamo un singolo punteggio di QI \(y_i\), la sua densità sotto una distribuzione normale con media \(\mu\) e deviazione standard \(\sigma=5\) è:

\[ p(y_i \mid \mu) = \text{Normal}(y_i \mid \mu, \sigma). \]

Quando abbiamo più osservazioni indipendenti, la plausibilità complessiva si ottiene moltiplicando le singole densità. La verosimiglianza congiunta diventa quindi:

\[ L(\mu) = \prod_{i=1}^n p(y_i \mid \mu, \sigma). \]

In altre parole, \(L(\mu)\) è alta se tutti i dati sono coerenti con un certo valore della media \(\mu\), mentre è bassa se i dati sono difficilmente spiegabili da quel valore.

Implementazione in R.

Per stimare la distribuzione a posteriori \(p(\mu \mid y)\) useremo il metodo della griglia. Definiamo prima una griglia di valori possibili per \(\mu\):

mu_griglia <- seq(120, 150, length.out = 100)
head(mu_griglia)
#> [1] 120.0 120.3 120.6 120.9 121.2 121.5

Poi calcoliamo la verosimiglianza congiunta per ciascun valore di \(\mu\). Per chiarezza usiamo un ciclo esplicito:

likelihood <- numeric(length(mu_griglia))  # vettore vuoto
for (j in seq_along(mu_griglia)) {
  mu <- mu_griglia[j]
  likelihood[j] <- prod(dnorm(campione, mean = mu, sd = 5))
}

Nell’istruzione:

likelihood[j] <- prod(dnorm(campione, mean = mu, sd = 5))

accade quanto segue:

  • dnorm(campione, mean = mu, sd = 5) calcola la densità normale per ciascuna delle 10 osservazioni del campione (ad es. 127, 129, 138, …).
  • prod(...) moltiplica tutte queste densità, costruendo la verosimiglianza congiunta per quel valore specifico di \(\mu\).
  • Questo procedimento viene ripetuto per ogni valore della griglia mu_griglia, così otteniamo una versione discretizzata della funzione di verosimiglianza.

Supponiamo di credere, prima di osservare i dati, che la media sia centrata su 140, con incertezza (deviazione standard 3). Formalmente: \(\mu \sim \mathcal{N}(140, 3^2)\).

prior <- dnorm(mu_griglia, mean = 140, sd = 3)

Il teorema di Bayes ci dice che la distribuzione a posteriori di \(\mu\), dato il campione osservato \(y\), è proporzionale al prodotto tra verosimiglianza e prior:

\[ p(\mu \mid y) \propto p(y \mid \mu) \cdot p(\mu). \]

In pratica, ogni valore di \(\mu\) riceve un “peso” che riflette sia la sua plausibilità secondo i dati (verosimiglianza), sia la sua plausibilità iniziale (prior).

Nel nostro approccio basato su griglia, otteniamo una versione discretizzata e non ancora normalizzata della distribuzione a posteriori:

posterior <- likelihood * prior

Per trasformarla in una vera distribuzione di probabilità discreta (cioè una distribuzione che somma a 1), normalizziamo dividendo per la somma totale di tutti i valori:

posterior <- posterior / sum(posterior)

Dopo questa operazione, posterior rappresenta una distribuzione di massa di probabilità sulla griglia di valori considerati per \(\mu\):

  • ogni elemento di posterior indica la probabilità a posteriori associata a un certo valore della griglia;
  • la somma di tutte queste probabilità è esattamente pari a 1.

Campionamento dalla posterior (Monte Carlo).

Una volta ottenuta la posterior discreta sulla griglia (mu_griglia) — cioè posterior con somma = 1 — possiamo stimare media e intervalli estraendo campioni proporzionali a tali probabilità.

# Sicurezza: la posterior deve sommare a 1
stopifnot(abs(sum(posterior) - 1) < 1e-10)

set.seed(123)
samples <- sample(mu_griglia, size = 10000, replace = TRUE, prob = posterior)

# Sommari Monte Carlo
mc_mean <- mean(samples)                       # media a posteriori (MC)
mc_ci94 <- quantile(samples, c(0.03, 0.97))    # intervallo di credibilità al 94%

mc_mean
#> [1] 132.6
mc_ci94
#>    3%   97% 
#> 130.0 135.2
  • La media a posteriori è la stima puntuale di \(\mu\) che integra dati e prior.
  • L’intervallo di credibilità al 94% (code uguali) contiene i valori di \(\mu\) che, a posteriori, coprono il 94% della probabilità.

Nota: essendo un metodo stocastico, i risultati variano leggermente a ogni estrazione (il seme stabilizza l’esempio). Con più campioni Monte Carlo, la variabilità si riduce.

Sommari a posteriori tramite integrazione su griglia (somme pesate).

In modo deterministico, possiamo calcolare media, varianza e quantili direttamente dalla griglia usando i pesi w = posterior (già normalizzati).

# Pesi e ascisse
w     <- posterior          # massa discreta sulla griglia (somma = 1)
theta <- mu_griglia

# Moda (MAP sulla griglia)
mode_hat <- theta[ which.max(w) ]

# Media e varianza (somme pesate)
mean_hat <- sum(theta * w)
var_hat  <- sum((theta - mean_hat)^2 * w)

# Intervallo al 94% per code uguali dalla CDF discreta
cdf  <- cumsum(w)
# Quantili come primi punti della griglia che raggiungono le probabilità target
q_lower <- theta[ which(cdf >= 0.03)[1] ]
q_upper <- theta[ which(cdf >= 0.97)[1] ]
ci94_det <- c(q_lower, q_upper)

data.frame(
  Stima  = c("Moda", "Media", "Varianza", "CI94_lower", "CI94_upper"),
  Valore = c(mode_hat, mean_hat, var_hat, ci94_det[1], ci94_det[2])
)
#>        Stima  Valore
#> 1       Moda 132.424
#> 2      Media 132.565
#> 3   Varianza   1.957
#> 4 CI94_lower 130.000
#> 5 CI94_upper 135.152
  • Poiché la griglia è equispaziata e posterior è già normalizzata, non serve alcun fattore \(\delta\).
  • L’uso della CDF discreta (cumsum) rende i quantili trasparenti senza bisogno di interpolazione.

Verifica dell’equivalenza (MC vs. griglia).

I due approcci devono dare risultati coerenti (a parte il rumore Monte Carlo):

abs(mean_hat - mc_mean)          # differenza media (attesa piccola)
#> [1] 5.27e-06
ci94_det
#> [1] 130.0 135.2
mc_ci94
#>    3%   97% 
#> 130.0 135.2

Visualizzazione.

Per capire come prior, verosimiglianza e posterior si combinano, è utile metterle sullo stesso grafico. Ricorda però che la verosimiglianza non è una distribuzione (non somma a 1): per il solo confronto visivo la normalizziamo (ad esempio dividendo per la somma).

# Normalizzazione SOLO per il confronto grafico
likelihood_std <- likelihood / sum(likelihood)

Questo non influisce in alcun modo sui calcoli dei sommari, che si basano unicamente sulla posterior.

dat <- tibble(
  mu = mu_griglia,
  Prior = prior / sum(prior),        # normalizzata solo per il grafico
  Verosimiglianza = likelihood_std,  # "finta densità" per confronto visivo
  Posteriori = posterior             # già massa discreta (somma = 1)
)

long_data <- dat |>
  pivot_longer(
    cols = c("Prior", "Verosimiglianza", "Posteriori"),
    names_to = "Distribuzione",
    values_to = "Densità"
  )

# Alcuni riferimenti utili sul grafico
mu_prior   <- 140
mu_sample  <- mean(campione)
mu_post    <- sum(mu_griglia * posterior)  # media a posteriori sulla griglia

ggplot(long_data, aes(x = mu, y = Densità, color = Distribuzione)) +
  geom_line(size = 1.2) +
  geom_vline(xintercept = mu_prior,  linetype = "dotted") +
  geom_vline(xintercept = mu_sample, linetype = "dashed") +
  geom_vline(xintercept = mu_post,   linetype = "dotdash") +
  annotate("text", x = mu_prior,  y = 0.06, label = "media prior",    vjust = 1.3, angle = 90, size = 3) +
  annotate("text", x = mu_sample, y = 0.06, label = "media campione", vjust = 1.3, angle = 90, size = 3) +
  annotate("text", x = mu_post,   y = 0.06, label = "media posterior",vjust = 1.3, angle = 90, size = 3) +
  labs(
    title = "Prior, verosimiglianza (normalizzata per confronto)\ne posterior",
    x = expression(mu),
    y = "Densità\n(normalizzata per confronto)",
    color = NULL
  ) +
  theme(legend.position = "bottom")

Interpretazione del grafico

  • Prior: rappresenta la convinzione iniziale (qui centrata su 140).
  • Verosimiglianza (normalizzata): mostra dove i dati sono più plausibili; tende a concentrarsi vicino alla media campionaria.
  • Posterior: è il compromesso tra prior e verosimiglianza. Con campioni più grandi e/o dati più informativi, la posterior si concentra maggiormente intorno ai valori suggeriti dai dati e diventa più stretta (minore incertezza).

Nota: la normalizzazione grafica rende confrontabili le forme delle curve (non le altezze assolute). I sommari numerici (media, varianza, intervalli) non usano questa normalizzazione: si basano solo su posterior.

Nota sulla stabilità numerica.

Poiché la verosimiglianza è un prodotto di molte probabilità molto piccole, il calcolo diretto può portare a valori quasi zero (problema di underflow). Per evitarlo, conviene lavorare con i logaritmi: il prodotto diventa una somma, più stabile numericamente.

log_likelihood <- numeric(length(mu_griglia))
for (j in seq_along(mu_griglia)) {
  mu <- mu_griglia[j]
  log_likelihood[j] <- sum(dnorm(campione, mean = mu, sd = 5, log = TRUE))
}

log_prior <- dnorm(mu_griglia, mean = 140, sd = 3, log = TRUE)

# log-posterior non normalizzata
log_posterior <- log_likelihood + log_prior

# stabilizzazione numerica (log-sum-exp trick)
log_posterior_stab <- log_posterior - max(log_posterior)

# posterior normalizzata su scala naturale
posterior_log <- exp(log_posterior_stab)
posterior_log <- posterior_log / sum(posterior_log)

Il risultato finale coincide con quello ottenuto senza logaritmi, ma in modo molto più sicuro dal punto di vista computazionale:

set.seed(123)
samples_log <- sample(mu_griglia, size = 10000, replace = TRUE, prob = posterior_log)
set.seed(123)
samples_dir <- sample(mu_griglia, size = 10000, replace = TRUE, prob = posterior)

c(mean_log = mean(samples_log), mean_dir = mean(samples_dir))
#> mean_log mean_dir 
#>    132.6    132.6
rbind(
  log = quantile(samples_log, c(0.03, 0.97)),
  dir = quantile(samples_dir, c(0.03, 0.97))
)
#>      3%   97%
#> log 130 135.2
#> dir 130 135.2

Perché i logaritmi sono importanti? (intuizione numerica)

Con 10 osservazioni l’underflow può non manifestarsi, ma i prodotti diventano numeri minuscoli:

range(likelihood)
#> [1] 1.576e-46 1.678e-13
formatC(min(likelihood), format = "e", digits = 2)  # notazione scientifica
#> [1] "1.58e-46"

Se aumentiamo il numero di osservazioni (ad es. replicando il campione per illustrazione), il prodotto può diventare 0 in aritmetica floating-point:

# Esempio illustrativo: campione "grande" replicando i dati (solo per mostrare l'effetto numerico)
campione_big <- rep(campione, 5)  # ~50 osservazioni

lik_big <- numeric(length(mu_griglia))
for (j in seq_along(mu_griglia)) {
  mu <- mu_griglia[j]
  lik_big[j] <- prod(dnorm(campione_big, mean = mu, sd = 5))
}

any(lik_big == 0)                 # può risultare TRUE (underflow)
#> [1] FALSE
formatC(min(lik_big), "e", 2)     # estremamente piccolo o 0
#> [1] "9.735208991935287351450858710642349463203201527141e-230"

Con i logaritmi questo problema sparisce, perché:

  • il prodotto diventa una somma: \(\sum \log p(y_i \mid \mu)\);
  • evitiamo l’overflow quando esponenziamo;
  • la normalizzazione finale su scala naturale restituisce una posterior ben definita.

Riflessioni conclusive

L’analisi presentata ha illustrato il processo di aggiornamento bayesiano sia per distribuzioni a priori discrete che continue, con particolare attenzione ai casi speciali che ammettono soluzioni analitiche. Questi casi, come la combinazione di priori Beta con verosimiglianze binomiali, offrono esempi particolarmente istruttivi del meccanismo bayesiano in azione.

Il metodo basato su griglia si è rivelato uno strumento didatticamente efficace per approssimare le distribuzioni a posteriori. La sua implementazione segue un percorso logico chiaro: si discretizza lo spazio dei parametri, si valutano prior e verosimiglianza in ciascun punto della griglia, e infine si normalizzano i risultati per ottenere una distribuzione di probabilità valida. Questo approccio mantiene una trasparenza concettuale che lo rende particolarmente adatto all’introduzione dei principi bayesiani.

Tuttavia, l’utilità pratica del metodo a griglia diminuisce rapidamente con l’aumentare della complessità del modello. La cosiddetta maledizione della dimensionalità rende proibitivo il costo computazionale quando si lavora con più di pochi parametri. Un modello con appena 10 parametri, utilizzando una griglia di 100 punti per ciascuno, richiederebbe la valutazione di un numero astronomico di combinazioni (\(10^{20}\)), dimostrando chiaramente i limiti di scalabilità di questo approccio.

Questa constatazione spiega perché nella pratica statistica avanzata si preferiscano tecniche più sofisticate come il campionamento MCMC. Tali metodi, pur essendo concettualmente più complessi, offrono la flessibilità e l’efficienza necessarie per affrontare problemi reali di moderata e alta dimensionalità. Il metodo a griglia rimane comunque prezioso come strumento didattico e per analisi preliminari in contesti semplici, rappresentando una porta d’accesso comprensibile al ricco panorama dell’inferenza bayesiana.

In uno studio sulla percezione delle emozioni, un partecipante osserva 10 fotografie di volti arrabbiati. Deve indicare se il volto esprime rabbia o no. Ogni risposta può essere corretta (1) o errata (0).

I dati osservati del partecipante sono:

1, 0, 1, 1, 1, 0, 0, 1, 1, 1

→ Totale: 7 successi su 10 prove → \(y = 7\), \(n = 10\).

Obiettivo: stimare la probabilità \(\theta\) che il partecipante riconosca correttamente un volto arrabbiato, tenendo conto sia dei dati osservati sia di conoscenze pregresse.

Prior Informativo.

Supponiamo di voler adottare un approccio cautamente pessimistico sulle capacità iniziali del partecipante, basandoci su studi precedenti che indicano un riconoscimento della rabbia non sempre accurato, ad esempio mediamente intorno al 40% con moderata incertezza.

Per rappresentare questa convinzione, scegliamo come distribuzione a priori una Beta(4, 6):

  • Media: \(\mu = \frac{4}{4+6} = 0.4\)
  • Varianza: \(\frac{4 \cdot 6}{(10)^2 \cdot 11} = 0.0218\)

Questa prior concentra la massa di probabilità su valori inferiori a 0.5, ma lascia spazio anche a livelli di competenza superiori.

Calcolo della Distribuzione a Posteriori con il Metodo Basato su Griglia.

1. Griglia di valori per \(\theta\).

theta <- seq(0, 1, length.out = 1000)
head(theta)
#> [1] 0.000000 0.001001 0.002002 0.003003 0.004004 0.005005
tail(theta)
#> [1] 0.995 0.996 0.997 0.998 0.999 1.000

2. Calcolo della distribuzione a priori Beta(4, 6).

prior <- dbeta(theta, shape1 = 4, shape2 = 6)
prior <- prior / sum(prior)  # normalizzazione

Visualizzazione:

ggplot(data.frame(theta, prior), aes(x = theta, y = prior)) +
  geom_line(linewidth = 1.2) +
  labs(
    title = "Distribuzione a Priori Informativa: Beta(4, 6)",
    x = expression(theta),
    y = "Densità (normalizzata)"
  ) +
  theme_minimal(base_size = 14)

3. Calcolo della verosimiglianza per 7 successi su 10.

likelihood <- dbinom(7, size = 10, prob = theta)
likelihood <- likelihood / sum(likelihood)  # normalizzazione

Visualizzazione:

ggplot(data.frame(theta, likelihood), aes(x = theta, y = likelihood)) +
  geom_line(linewidth = 1.2) +
  labs(
    title = "Funzione di Verosimiglianza: 7 successi su 10",
    x = expression(theta),
    y = "Densità (normalizzata)"
  ) +
  theme_minimal(base_size = 14)

4. Calcolo della distribuzione a posteriori.

unnormalized_posterior <- prior * likelihood
posterior <- unnormalized_posterior / sum(unnormalized_posterior)

5. Visualizzazione congiunta: prior, likelihood e posteriori.

data <- data.frame(theta, prior, likelihood, posterior)

# Imposta i livelli desiderati con nomi leggibili
long_data <- pivot_longer(
  data,
  cols = c("prior", "likelihood", "posterior"),
  names_to = "distribution",
  values_to = "density"
) |>
  mutate(distribution = factor(
    distribution,
    levels = c("prior", "likelihood", "posterior"),
    labels = c("A Priori", "Verosimiglianza", "A Posteriori")
    )
  )

ggplot(
  long_data, 
  aes(x = theta, y = density, color = distribution)
  ) +
  geom_line(linewidth = 1.2) +
  labs(
    title = "Distribuzioni Bayesiane: A Priori, Verosimiglianza e A Posteriori",
    x = expression(theta),
    y = "Densità (normalizzata)",
    color = NULL
  ) +
  theme(legend.position = "bottom")

Riepilogo:

  • la prior (Beta(4,6)) riflette una convinzione iniziale più scettica;
  • la verosimiglianza è centrata su \(\theta = 0.7\), corrispondente a 7 successi su 10;
  • la posteriori media tra prior e dati, ma si sposta chiaramente verso destra, evidenziando l’effetto aggiornamento bayesiano.

Questo esempio mostra come l’approccio bayesiano:

  • integra in modo trasparente dati individuali e credenze pregresse;
  • produce una stima personalizzata della capacità del partecipante;
  • permette di quantificare l’incertezza in modo completo, tramite la distribuzione a posteriori.

Quantità a Posteriori.

Media:

posterior_mean <- sum(theta * posterior)
posterior_mean
#> [1] 0.55

Deviazione standard:

posterior_sd <- sqrt(sum((theta^2) * posterior) - posterior_mean^2)
posterior_sd
#> [1] 0.1086

Moda:

posterior_mode <- theta[which.max(posterior)]
posterior_mode
#> [1] 0.5556

Intervallo di credibilità al 94%:

samples <- sample(theta, size = 10000, replace = TRUE, prob = posterior)
quantile(samples, probs = c(0.03, 0.97))
#>     3%    97% 
#> 0.3433 0.7528

Questo esercizio mostra come:

  • l’informazione pregressa può essere incorporata in modo trasparente in un modello bayesiano;
  • la posteriori riflette una combinazione tra dati osservati e conoscenze precedenti.

In uno studio sull’analisi delle pratiche di trasparenza e riproducibilità nella ricerca in psicologia, Hardwicke et al. (2022) hanno riportato che la condivisione dei materiali di ricerca è stata rilevata nel 14% dei casi (26 su 183 studi), con un intervallo di confidenza al 95% pari a [10%, 19%]. Questo suggerisce che la condivisione di materiali è rara.

Ispirandoti ai risultati di questo studio, costruisci una distribuzione a priori per la probabilità \(\theta\) che uno studio condivida i materiali di ricerca. Per semplicità, discretizza \(\theta\) in 10 livelli equispaziati: \(0.05, 0.15, 0.25, 0.35, 0.45, 0.55, 0.65, 0.75, 0.85, 0.95\).

Attribuisci le seguenti probabilità a priori ai 10 livelli, basandoti sull’informazione che la condivisione dei materiali è un evento raro ma non trascurabile: \(0.05, 0.20, 0.30, 0.15, 0.10, 0.08, 0.05, 0.03, 0.02, 0.02\).

Supponiamo che siano stati osservati 20 studi su 100 che hanno condiviso i materiali di ricerca. Calcola la distribuzione a posteriori utilizzando il metodo basato su griglia. Calcola la media della distribuzione a posteriori e l’intervallo di credibilità al 89%.

# Definizione dei possibili valori di theta (probabilità discreta)
theta <- seq(0.05, 0.95, by = 0.10)

# Definizione della distribuzione a priori
prior <- c(0.05, 0.20, 0.30, 0.15, 0.10, 0.08, 0.05, 0.03, 0.02, 0.02)

# Normalizzazione della prior (se necessario, ma in questo caso già normalizzata)
prior <- prior / sum(prior)

# Dati osservati
successi <- 20  # studi che hanno condiviso materiali
n <- 100        # studi totali

# Calcolo della verosimiglianza usando la distribuzione binomiale
likelihood <- dbinom(successi, size = n, prob = theta)

# Calcolo della distribuzione a posteriori (applicazione del teorema di Bayes)
posterior <- likelihood * prior
posterior <- posterior / sum(posterior)  # Normalizzazione

# Calcolo della media della distribuzione a posteriori
posterior_mean <- sum(theta * posterior)

# Calcolo dell'intervallo di credibilità al 89%
cdf <- cumsum(posterior)  # Distribuzione cumulativa
lower_bound <- theta[which.min(abs(cdf - 0.055))]  # 5.5% quantile
upper_bound <- theta[which.min(abs(cdf - 0.945))]  # 94.5% quantile

# Output dei risultati
list(
  posterior_mean = posterior_mean,
  credibility_interval_89 = c(lower_bound, upper_bound)
)

L’obiettivo di questo esercizio è applicare il metodo basato su griglia per stimare la distribuzione a posteriori di una proporzione, utilizzando dati dalla Scala della Rete Sociale di Lubben (LSNS-6). Si assume che un punteggio LSNS-6 superiore a una soglia prefissata indichi isolamento sociale. Il compito è:

  1. Scegliere una soglia per classificare i partecipanti in due gruppi (isolati vs. non isolati), garantendo che la proporzione osservata non sia inferiore a 0.1 o superiore a 0.9.
  2. Calcolare la distribuzione a posteriori della proporzione usando un’approssimazione discreta su una griglia di valori.
  3. Determinare l’intervallo di credibilità all’89%.
  4. Interpretare i risultati.

Consegna: caricare su Moodle il file .qmd compilato in pdf.

Dati e Modellizzazione

Si assume che i dati siano rappresentati da una variabile binaria \(y\), con \(y = 1\) per individui classificati come isolati e \(y = 0\) altrimenti. Supponiamo che su un campione di \(n\) individui, \(s\) siano isolati.

Definiamo il modello statistico:

\[ y_i \sim \text{Bernoulli}(\theta) \]

con:

  • \(y_i \in \{0,1\}\) per \(i=1,\dots,n\),
  • \(\theta\) proporzione di individui isolati nella popolazione.

La distribuzione a priori su \(\theta\) è scelta come \(\text{Beta}(2,2)\), che rappresenta una conoscenza iniziale moderata e non estrema.

Metodo basato su griglia

Il metodo a griglia approssima la distribuzione a posteriori calcolando la probabilità per una serie di valori discreti di \(\theta\).

  1. Definire una griglia di valori per \(\theta\):

    theta <- seq(0, 1, length.out = 100)

  2. Calcolare la distribuzione a priori:

    prior <- dbeta(theta, 2, 2)
prior <- prior / sum(prior)  # Normalizzazione

  3. Calcolare la verosimiglianza:

    likelihood <- dbinom(s, size = n, prob = theta)
likelihood <- likelihood / sum(likelihood)  # Normalizzazione

  4. Calcolare la distribuzione a posteriori:

    posterior <- prior * likelihood
posterior <- posterior / sum(posterior)  # Normalizzazione

** Calcolo dell’intervallo di credibilità all’89%**

L’intervallo di credibilità è calcolato come l’intervallo che contiene il 89% della probabilità a posteriori.

ci_89 <- quantile(sample(theta, size = 10000, prob = posterior, replace = TRUE), probs = c(0.055, 0.945))
ci_89

Interpretazione dei risultati

  1. Valore atteso e moda a posteriori:

    mean_theta <- sum(theta * posterior)
mode_theta <- theta[which.max(posterior)]
    • Il valore atteso fornisce una stima puntuale di \(\theta\).
    • La moda indica il valore più probabile della proporzione di isolamento sociale.

  2. Intervallo di credibilità:

    • L’89% della probabilità a posteriori cade tra i valori dell’intervallo di credibilità.
    • Se l’intervallo è stretto, c’è maggiore certezza sulla proporzione stimata.
    • Se l’intervallo è ampio, vi è maggiore incertezza sulla proporzione.
sessionInfo()
#> R version 4.5.1 (2025-06-13)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS Sequoia 15.6.1
#> 
#> 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] HDInterval_0.2.4      pillar_1.11.0         tinytable_0.11.0     
#>  [4] patchwork_1.3.1       ggdist_3.3.3          tidybayes_3.0.7      
#>  [7] bayesplot_1.13.0      ggplot2_3.5.2         reliabilitydiag_0.2.1
#> [10] priorsense_1.1.0      posterior_1.6.1       loo_2.8.0            
#> [13] rstan_2.32.7          StanHeaders_2.32.10   brms_2.22.0          
#> [16] Rcpp_1.1.0            sessioninfo_1.2.3     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] svUnit_1.0.6         tidyselect_1.2.1     farver_2.1.2        
#>  [4] fastmap_1.2.0        TH.data_1.1-3        tensorA_0.36.2.1    
#>  [7] pacman_0.5.1         digest_0.6.37        timechange_0.3.0    
#> [10] estimability_1.5.1   lifecycle_1.0.4      survival_3.8-3      
#> [13] magrittr_2.0.3       compiler_4.5.1       rlang_1.1.6         
#> [16] tools_4.5.1          yaml_2.3.10          knitr_1.50          
#> [19] labeling_0.4.3       bridgesampling_1.1-2 htmlwidgets_1.6.4   
#> [22] curl_6.4.0           pkgbuild_1.4.8       RColorBrewer_1.1-3  
#> [25] abind_1.4-8          multcomp_1.4-28      withr_3.0.2         
#> [28] purrr_1.1.0          grid_4.5.1           stats4_4.5.1        
#> [31] colorspace_2.1-1     xtable_1.8-4         inline_0.3.21       
#> [34] emmeans_1.11.2       scales_1.4.0         MASS_7.3-65         
#> [37] cli_3.6.5            mvtnorm_1.3-3        rmarkdown_2.29      
#> [40] ragg_1.4.0           generics_0.1.4       RcppParallel_5.1.10 
#> [43] cachem_1.1.0         stringr_1.5.1        splines_4.5.1       
#> [46] parallel_4.5.1       vctrs_0.6.5          V8_6.0.5            
#> [49] Matrix_1.7-3         sandwich_3.1-1       jsonlite_2.0.0      
#> [52] arrayhelpers_1.1-0   systemfonts_1.2.3    glue_1.8.0          
#> [55] codetools_0.2-20     distributional_0.5.0 lubridate_1.9.4     
#> [58] stringi_1.8.7        gtable_0.3.6         QuickJSR_1.8.0      
#> [61] htmltools_0.5.8.1    Brobdingnag_1.2-9    R6_2.6.1            
#> [64] textshaping_1.0.1    rprojroot_2.1.0      evaluate_1.0.4      
#> [67] lattice_0.22-7       backports_1.5.0      memoise_2.0.1       
#> [70] broom_1.0.9          snakecase_0.11.1     rstantools_2.4.0    
#> [73] coda_0.19-4.1        gridExtra_2.3        nlme_3.1-168        
#> [76] checkmate_2.3.2      xfun_0.52            zoo_1.8-14          
#> [79] pkgconfig_2.0.3

Bibliografia

Albert, J., & Hu, J. (2019). Probability and Bayesian Modeling. CRC Press.
Hardwicke, T. E., Thibault, R. T., Kosie, J. E., Wallach, J. D., Kidwell, M. C., & Ioannidis, J. P. (2022). Estimating the prevalence of transparency and reproducibility-related research practices in psychology (2014–2017). Perspectives on Psychological Science, 17(1), 239–251.