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

In questo capitolo imparerai a

• 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.

Prerequisiti

• Leggere il settimo capitolo del testo di Albert & Hu (2019).

Preparazione del Notebook

here::here("code", "_common.R") |> 
  source()

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

47.1 Introduzione

L’inferenza bayesiana offre un quadro formale per aggiornare le nostre convinzioni riguardo a un parametro incognito, come la probabilità \(\theta\) di un evento, alla luce di nuove osservazioni. Questo approccio è basato sul teorema di Bayes, che combina in modo coerente le conoscenze iniziali—espresse attraverso una distribuzione a priori—con l’evidenza fornita dai dati, sintetizzata dalla funzione di verosimiglianza. Il risultato è una distribuzione a posteriori, che rappresenta le nostre credenze rivedute dopo aver considerato sia l’ipotesi iniziale sia i dati osservati.

L’inferenza bayesiana si distingue per la sua flessibilità nel trattare l’incertezza: invece di fornire stime puntuali, lavora con distribuzioni di probabilità, permettendo di quantificare in modo esplicito la fiducia associata a ciascun valore possibile del parametro. Tale caratteristica la rende particolarmente utile in contesti in cui le informazioni sono incomplete o soggette a revisione, come spesso accade nella ricerca psicologica.

In questo capitolo esploreremo l’inferenza bayesiana attraverso il modello binomiale, particolarmente adatto per analizzare dati dicotomici. Inizieremo considerando un semplice caso con distribuzione a priori discreta, dove il parametro \(\theta\) (probabilità di successo) può assumere solo valori predefiniti. Questo approccio ci permetterà di osservare come le nostre credenze su \(\theta\) si modifichino gradualmente con l’arrivo di nuovi dati. Successivamente, estenderemo il discorso al caso più generale delle distribuzioni a priori continue, che meglio si adattano alla modellizzazione di problemi reali. L’obiettivo è fornire una comprensione intuitiva ma rigorosa del processo dell’aggiornamento bayesiano, evidenziandone sia la logica sottostante sia le potenzialità applicative.

47.2 Verosimiglianza Binomiale

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

\[ 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).

47.3 Esempio: Riconoscimento di Parole Emozionanti

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.

47.3.1 Obiettivo inferenziale

Vogliamo stimare il parametro \(\theta\), che rappresenta la probabilità che un adulto, estratto a caso dalla popolazione, riconosca correttamente parole emotive. In un contesto bayesiano, il modello binomiale ci permette di stimare \(\theta\) a partire dai dati osservati e di quantificare l’incertezza associata a questa stima. In particolare, la verosimiglianza dei dati osservati (22 successi su 30 prove) può essere espressa attraverso la funzione di massa di probabilità binomiale, che assegna una probabilità a ogni possibile valore di \(\theta\) compreso tra 0 e 1.

L’analisi bayesiana di questo problema prevede la combinazione di queste informazioni con una distribuzione a priori che rappresenta le nostre credenze iniziali sull’abilità inibitoria del partecipante. Questo approccio ci consentirà di ottenere una distribuzione a posteriori per \(\theta\) che sintetizza sia le informazioni a priori sia l’evidenza fornita dai dati osservati.

47.4 Metodo Basato su Griglia per l’Aggiornamento Bayesiano

Il metodo basato su griglia rappresenta un approccio computazionalmente semplice per ottenere la distribuzione a posteriori quando non è disponibile una soluzione analitica diretta. Questo metodo risulta particolarmente utile per illustrare concretamente il processo di aggiornamento bayesiano.

Il procedimento inizia definendo un intervallo plausibile per il parametro \(\theta\), tipicamente compreso tra 0 e 1 nel caso di probabilità. All’interno di questo intervallo, viene costruita una griglia costituita da un insieme finito di valori equidistanti (ad esempio, \(\theta\) = 0, 0.01, 0.02, …, 1). Per ciascun valore \(\theta\) nella griglia, si calcola il prodotto della verosimiglianza dei dati osservati (data \(\theta\)) e della densità a priori corrispondente.

I valori così ottenuti vengono poi normalizzati dividendo ciascuno di essi per la somma totale dei prodotti calcolati su tutta la griglia. Questa operazione garantisce che l’area sotto la curva risultante sia pari a 1, trasformando così i valori in una vera distribuzione di probabilità. La distribuzione risultante rappresenta un’approssimazione discreta della vera distribuzione a posteriori continua.

Sebbene approssimato, questo metodo offre diversi vantaggi didattici: mantiene trasparente ogni fase del calcolo, permette di visualizzare direttamente l’effetto dell’aggiornamento bayesiano e non richiede conoscenze avanzate di metodi computazionali. La precisione dell’approssimazione può essere migliorata aumentando la densità dei punti nella griglia, a scapito però di un maggior costo computazionale.

Nel caso del nostro esempio con 22 successi su 30 prove, applicheremo questo metodo utilizzando come a priori una distribuzione uniforme, dimostrando come l’osservazione dei dati modifichi la nostra comprensione della probabilità \(\theta\) di riconoscimento corretto di parole con contenuto emotivo positivo.

47.5 Aggiornamento Bayesiano con Priori Discreta

47.5.1 Costruzione della Priori

Nel caso del compito di compito di memoria a lungo termine considerato, il parametro \(theta\) rappresenta la probabilità di un riconoscimento corretto di parole con contenuto emotivo positivo. In assenza di informazioni preliminari specifiche, potremmo inizialmente considerare tutti i valori di \(\theta\) 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.

47.5.2 Aggiornamento con i Dati

Supponiamo di osservare 22 inibizioni corrette su 30 prove (\(y\)=22, \(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 maggiore probabilità a posteriori.

47.5.3 Interpretazione

Questo processo dimostra come l’osservazione dei dati modifichi le nostre credenze iniziali. La distribuzione a posteriori si sposta gradualmente dalla configurazione iniziale uniforme verso una forma più informata, concentrandosi attorno ai valori di \(\theta\) più compatibili con i dati osservati.

In particolare, dopo aver osservato 22 successi su 30 prove, valori come \(\theta\) = 0.7 o 0.75 diventano più plausibili rispetto a valori estremi come \(\theta\) = 0.2 o 0.9. Questo riflette il modo in cui l’evidenza empirica aggiorna le nostre aspettative iniziali.

L’uso di una griglia discreta fornisce un’approssimazione della vera distribuzione continua, dove la precisione dell’approssimazione migliora quanto più fitta è la griglia di punti considerata. Questo approccio permette di visualizzare in modo concreto il meccanismo di aggiornamento bayesiano, mostrando come a ogni nuova osservazione corrisponda una revisione della nostra conoscenza del parametro \(\theta\).

47.5.4 Implementazione in R

47.5.4.1 Definizione dei parametri

Iniziamo creando una griglia discreta per il parametro θ, che rappresenta la probabilità di successo nel nostro compito Go/No-Go:

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

47.5.4.2 Distribuzione a priori uniforme

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à")

47.5.4.3 Distribuzione a priori informativa

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à")

47.5.4.4 Calcolo della verosimiglianza

Per i nostri dati (22 successi su 30 prove):

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

47.5.4.5 Calcolo della distribuzione a posteriori

Combiniamo priori e verosimiglianza:

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

47.5.4.6 Statistiche descrittive

Calcoliamo alcune quantità riassuntive:

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.

La metodologia si distingue per la sua flessibilità applicativa, permettendo di adeguarsi a diverse esigenze analitiche. È possibile modificare facilmente i parametri osservati, come il numero di successi e prove registrate, per adattare l’analisi a diversi set di dati. La distribuzione a priori può essere calibrata in base alle specifiche ipotesi teoriche del ricercatore, offrendo la possibilità di incorporare conoscenze pregresse di varia natura. Un ulteriore vantaggio operativo consiste nella possibilità di aumentare la precisione delle stime semplicemente aggiungendo punti alla griglia di valori considerati per \(\theta\), ottenendo così un’approssimazione più fedele della distribuzione continua. Questa combinazione di adattabilità e precisione rende l’approccio particolarmente efficace sia nelle fasi iniziali di analisi esplorativa, dove è necessario testare rapidamente diverse configurazioni, sia in contesti sperimentali controllati che richiedono modelli più sofisticati.

47.6 Aggiornamento Bayesiano con una Distribuzione a Priori Continua

Passiamo ora a un’estensione naturale del metodo discreto: l’utilizzo di una distribuzione a priori continua. Questo approccio permette di modellare la probabilità di successo \(\theta\) come una variabile continua definita sull’intervallo \([0, 1]\), consentendo maggiore precisione e flessibilità. In questo contesto, la distribuzione Beta rappresenta una scelta ideale per la prior su una proporzione, grazie alla sua versatilità e alla sua coniugatezza con la distribuzione binomiale.

47.6.1 Scelta della Distribuzione a Priori

Iniziamo con una distribuzione Beta simmetrica, \(\text{Beta}(2, 2)\), che riflette una credenza iniziale moderatamente incerta, ma senza una preferenza marcata per valori bassi o alti di \(\theta\). Successivamente, considereremo anche un caso più informativo: \(\text{Beta}(2, 5)\), che assegna maggiore probabilità a valori di \(\theta\) vicini a zero, riflettendo un’aspettativa più pessimista.

47.6.2 Costruzione della Densità a Priori

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

Visualizziamo la forma della distribuzione:

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à")

47.6.3 Verosimiglianza per un Esperimento Binomiale

Supponiamo che un partecipante abbia ottenuto 6 successi su 9 prove in un compito dicotomico (es. compito No-Go). La verosimiglianza associata a ciascun valore di \(\theta\) è calcolata come:

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

47.6.4 Calcolo della Distribuzione a Posteriori

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

47.6.5 Distribuzione a Priori Non Simmetrica

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

47.6.6 Distribuzione a Posteriori

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

47.6.7 Quantità a Posteriori

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

47.6.8 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.

47.7 Estensione: Il Metodo Basato su Griglia per il Modello con Verosimiglianza Gaussiana

Finora abbiamo applicato il metodo basato su griglia al caso della stima di una proporzione in un modello beta-binomiale. Tuttavia, lo stesso principio si estende in modo naturale anche a modelli con parametri continui e verosimiglianza diversa dalla binomiale, come nel caso della distribuzione normale.

In particolare, 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.

L’uso della griglia permette di approssimare la distribuzione a posteriori in modo trasparente e intuitivo, anche quando il modello non ha una soluzione analitica semplice.

47.7.1 Esempio: Stima Bayesiana della Media del QI

Immaginiamo di condurre uno studio su bambini ad alto potenziale cognitivo. Per semplicità, ipotizziamo che il QI in questa popolazione segua una distribuzione normale con deviazione standard nota pari a 5, mentre la media \(\mu\) è incognita e rappresenta il parametro di interesse.

Supponiamo di osservare i seguenti 10 punteggi di QI (simulati da una normale con media 130 e sd = 5):

set.seed(123)
campione <- round(rnorm(10, mean = 130, sd = 5))

47.7.2 Definizione della Griglia

Per applicare il metodo bayesiano su griglia, costruiamo una sequenza di valori plausibili per \(\mu\), ad esempio:

mu_griglia <- seq(110, 150, length.out = 100)

47.7.3 Calcolo della Verosimiglianza

La verosimiglianza per ciascun valore di \(\mu\) è il prodotto delle densità normali dei dati osservati, assumendo \(\sigma = 5\):

likelihood <- sapply(mu_griglia, function(mu) {
  prod(dnorm(campione, mean = mu, sd = 5))
})

47.7.4 Prior Informativa

Supponiamo ora di avere una convinzione iniziale secondo cui \(\mu\) è centrata intorno a 140, con una deviazione standard pari a 3. Questa informazione può essere rappresentata con una distribuzione a priori gaussiana \(\mathcal{N}(140, 3^2)\):

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

Combinando questa prior con la verosimiglianza calcolata in precedenza, otteniamo la distribuzione a posteriori:

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

Questa posteriori riflette un compromesso tra i dati osservati e la credenza iniziale. Se i dati sono informativi (campione sufficientemente grande), l’influenza della prior si riduce; viceversa, con pochi dati, la prior può giocare un ruolo più rilevante.

47.7.5 Campionamento dalla Posteriori

Possiamo generare un campione dalla distribuzione a posteriori discreta per calcolare quantità riassuntive e costruire intervalli di credibilità:

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

Quantità riassuntive:

mean(samples)                         # media a posteriori
#> [1] 132.6
quantile(samples, c(0.03, 0.97))      # intervallo di credibilità al 94%
#>    3%   97% 
#> 129.8 135.1

47.7.6 Visualizzazione congiunta: Prior, Verosimiglianza, Posteriori

Per confrontare le tre distribuzioni, normalizziamo anche la verosimiglianza:

likelihood_std <- likelihood / sum(likelihood)

dat <- tibble(
  mu = mu_griglia,
  Prior = prior / sum(prior),
  Verosimiglianza = likelihood_std,
  Posteriori = posterior
)

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

Grafico:

ggplot(long_data, aes(x = mu, y = Densità, color = Distribuzione)) +
  geom_line(size = 1.2) +
  labs(
    title = "Confronto tra Prior, Verosimiglianza e Posteriori",
    x = expression(mu),
    y = "Densità (normalizzata)",
    color = NULL
  ) +
  theme(legend.position = "bottom")

47.7.7 Interpretazione

• La prior riflette una credenza iniziale centrata su \(\mu = 140\), coerente con l’ipotesi che i bambini siano plusdotati.
• La verosimiglianza riflette i dati osservati, tendendo a concentrare la densità intorno alla media campionaria.
• La posteriori media tra queste due fonti di informazione, spostandosi più o meno verso la prior a seconda del numero di osservazioni e della forza dei dati.

Questo esempio mostra chiaramente come l’inferenza bayesiana combini in modo coerente dati e credenze preesistenti, rendendo il processo inferenziale trasparente e aggiornabile.

47.7.8 Nota su Stabilità Numerica: Uso dei Logaritmi

Per evitare problemi di underflow, è possibile calcolare i logaritmi delle densità:

log_likelihood <- sapply(mu_griglia, function(mu) {
  sum(dnorm(campione, mean = mu, sd = 5, log = TRUE))
})
log_prior <- dnorm(mu_griglia, mean = 140, sd = 3, log = TRUE)

log_posterior <- log_likelihood + log_prior
log_posterior <- log_posterior - max(log_posterior)  # stabilizzazione
posterior <- exp(log_posterior) / sum(exp(log_posterior))

Le quantità riassuntive che otteniamo in questo modo coincidono con quelle trovate in precedenza:

set.seed(123)
samples <- sample(mu_griglia, size = 10000, replace = TRUE, prob = posterior)
mean(samples)                         # media a posteriori
#> [1] 132.6
quantile(samples, c(0.03, 0.97))      # intervallo di credibilità al 94%
#>    3%   97% 
#> 129.8 135.1

47.7.9 Estensione: Parametri Multipli

Nel caso in cui anche \(\sigma\) sia ignota, possiamo estendere la griglia a due dimensioni:

• \(\mu \in [110, 150]\)
• \(\sigma \in [1, 10]\)

e calcolare la verosimiglianza e la prior su ciascuna coppia \((\mu, \sigma)\). Tuttavia, il numero di combinazioni cresce rapidamente, evidenziando la maledizione della dimensionalità. In tali casi, è consigliabile usare metodi di approssimazione come l’MCMC.

47.8 Riflessioni conclusive

In questo capitolo abbiamo esaminato l’aggiornamento bayesiano nel caso di distribuzioni a priori discrete, con un breve accenno alle prior continue. Mentre queste ultime richiedono solitamente metodi numerici per il calcolo della distribuzione a posteriori, esistono situazioni particolari—come l’inferenza su proporzioni con prior Beta e verosimiglianza binomiale—in cui la posterior è derivabile analiticamente. Questi casi saranno approfonditi nei capitoli successivi.

L’approccio utilizzato qui si basa sul metodo basato su griglia, una tecnica numerica che approssima la posterior discretizzando lo spazio dei parametri. Il procedimento è semplice: si definisce una griglia di valori plausibili per il parametro, si valutano prior e verosimiglianza in ciascun punto, si calcola la posterior non normalizzata e infine la si normalizza per ottenere una distribuzione di probabilità valida. Questo metodo fornisce risultati esatti (a meno dell’approssimazione della griglia) e non richiede strumenti computazionali avanzati.

Il metodo basato su griglia si applica anche a modelli continui con verosimiglianza gaussiana, come nel caso della stima della media. L’idea centrale resta invariata: valutare, per ciascun punto della griglia, il prodotto tra prior e verosimiglianza, e normalizzare. L’approccio mantiene una trasparenza didattica elevata, rendendolo ideale per introdurre l’inferenza bayesiana in contesti psicologici anche con variabili continue.

Tuttavia, l’utilità del metodo basato su griglia è limitata dalla maledizione della dimensionalità: mentre per uno o due parametri il metodo rimane efficiente, all’aumentare del loro numero il costo computazionale diventa rapidamente proibitivo. Basti pensare che, con una griglia di 100 punti per parametro, un modello a 10 parametri richiederebbe la valutazione di \(10^{20}\) combinazioni, rendendo l’approccio impraticabile. Per questo, nella statistica bayesiana applicata—dove i modelli spesso coinvolgono decine o centinaia di parametri—si preferiscono tecniche più efficienti, come il campionamento MCMC, che affronteremo in seguito.

In sintesi, il metodo a griglia rappresenta un’introduzione accessibile all’inferenza bayesiana, utile per comprendere i concetti fondamentali e applicabile a problemi semplici. Tuttavia, la sua scalabilità limitata ne riduce l’uso pratico, spingendo verso approcci più sofisticati quando si lavora con modelli realistici.

Esercizi

Esercizio

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, color = okabe_ito_palette[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, color = okabe_ito_palette[1]) +
  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
  ) +
  scale_color_manual(
    values = c(
      "A Priori" = okabe_ito_palette[2],       # blu
      "Verosimiglianza" = okabe_ito_palette[1],# arancione
      "A Posteriori" = okabe_ito_palette[3]    # verde
    )
  ) +
  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.

Problemi 1

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%.

Soluzioni 1

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

Problemi 2

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.

Soluzioni 2

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.

Informazioni sull’Ambiente di Sviluppo

sessionInfo()
#> R version 4.5.0 (2025-04-11)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS Sequoia 15.5
#> 
#> 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/Rome
#> tzcode source: internal
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#>  [1] HDInterval_0.2.4 thematic_0.1.7   MetBrewer_0.2.0  ggokabeito_0.1.0
#>  [5] see_0.11.0       gridExtra_2.3    patchwork_1.3.0  bayesplot_1.13.0
#>  [9] psych_2.5.3      scales_1.4.0     markdown_2.0     knitr_1.50      
#> [13] lubridate_1.9.4  forcats_1.0.0    stringr_1.5.1    dplyr_1.1.4     
#> [17] purrr_1.0.4      readr_2.1.5      tidyr_1.3.1      tibble_3.3.0    
#> [21] ggplot2_3.5.2    tidyverse_2.0.0  rio_1.2.3        here_1.0.1      
#> 
#> loaded via a namespace (and not attached):
#>  [1] generics_0.1.4     stringi_1.8.7      lattice_0.22-7    
#>  [4] hms_1.1.3          digest_0.6.37      magrittr_2.0.3    
#>  [7] evaluate_1.0.4     grid_4.5.0         timechange_0.3.0  
#> [10] RColorBrewer_1.1-3 fastmap_1.2.0      rprojroot_2.0.4   
#> [13] jsonlite_2.0.0     mnormt_2.1.1       cli_3.6.5         
#> [16] rlang_1.1.6        withr_3.0.2        tools_4.5.0       
#> [19] parallel_4.5.0     tzdb_0.5.0         pacman_0.5.1      
#> [22] vctrs_0.6.5        R6_2.6.1           lifecycle_1.0.4   
#> [25] htmlwidgets_1.6.4  pkgconfig_2.0.3    pillar_1.10.2     
#> [28] gtable_0.3.6       glue_1.8.0         xfun_0.52         
#> [31] tidyselect_1.2.1   rstudioapi_0.17.1  farver_2.1.2      
#> [34] htmltools_0.5.8.1  nlme_3.1-168       labeling_0.4.3    
#> [37] rmarkdown_2.29     compiler_4.5.0

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.