here::here("code", "_common.R") |>
source()
# Load packages
if (!requireNamespace("pacman")) install.packages("pacman")
pacman::p_load(mice)52 Modello coniugato Gamma-Poisson
52.1 Introduzione
In psicologia, le variabili di conteggio (\(y\)), che indicano il numero di occorrenze di un evento, trovano ampio impiego in diversi ambiti. Ad esempio, sono usate per quantificare la frequenza dei sintomi di un disturbo o per analizzare le frequenze delle parole negli studi di psicolinguistica. Queste variabili, assumendo valori discreti, richiedono modelli statistici specifici.
Questo capitolo si focalizza sulla stima del tasso medio di incidenza (\(\lambda_i\)) di tali eventi, ovvero sul numero medio di occorrenze per unità di misura. Adotteremo un approccio bayesiano, utilizzando il modello di Poisson per descrivere la distribuzione di probabilità delle variabili di conteggio. In particolare, esploreremo la derivazione analitica della distribuzione a posteriori del parametro \(\lambda_i\), considerando una distribuzione a priori Gamma. Successivamente, verificheremo la validità dei risultati analitici mediante simulazioni Monte Carlo.
52.2 Distribuzione di Poisson
La distribuzione di Poisson è un modello probabilistico utilizzato per descrivere il numero di eventi che si verificano in un intervallo di tempo o spazio fisso, partendo dall’assunto che tali eventi si verifichino con una frequenza media costante e in modo indipendente rispetto al tempo trascorso dall’ultimo evento. Se un dato \(y\) segue una distribuzione di Poisson con parametro \(\lambda\), allora la probabilità di osservare un singolo valore \(y_i\) è data da:
\[ f(y_i \mid \lambda) = \frac{e^{-\lambda} \lambda^{y_i}}{y_i!}, \]
dove \(\lambda > 0\) rappresenta la frequenza media di occorrenza degli eventi e \(y_i\) è il numero di eventi osservati. La distribuzione di Poisson ha la caratteristica che sia il valore atteso che la varianza di una variabile casuale \(Y\) che segue questa distribuzione sono pari a \(\lambda\), cioè \(E(Y) = \lambda\) e \(\text{Var}(Y) = \lambda\).
52.2.1 Simulazione
Per capire meglio come funziona la distribuzione di Poisson, immaginiamo un paziente con un disturbo ossessivo-compulsivo. Supponiamo che in media questo paziente ripeta un’azione compulsiva 2 volte ogni ora. In questo caso, il parametro della distribuzione di Poisson è λ = 2.
La probabilità di osservare esattamente \(k\) eventi in un’ora è calcolata dalla formula:
\[ f(k \mid \lambda) = \frac{e^{-\lambda} \lambda^k}{k!}. \]
Nel caso specifico con \(\lambda = 2\), le probabilità per i primi valori di \(k\) sono:
- La probabilità di osservare 0 eventi in un’ora è \(\frac{e^{-2} \cdot 2^0}{0!} = e^{-2} \approx 0{.}1353\).
- La probabilità di osservare 1 evento in un’ora è \(\frac{e^{-2} \cdot 2^1}{1!} = 2e^{-2} \approx 0{.}2707\).
- La probabilità di osservare 2 eventi in un’ora è \(\frac{e^{-2} \cdot 2^2}{2!} = 2e^{-2} \approx 0{.}2707\).
- E così via per \(k = 3\), \(k = 4\), \(\dots\)
Questo esempio illustra come la distribuzione di Poisson possa essere utilizzata per modellare il numero di eventi rari che si verificano in un intervallo temporale fisso, con una frequenza media nota.
Svolgiamo ora i calcoli utilizzando la funzione dpois in R:
lam_true <- 2
# Creazione di un vettore di valori da 0 a 9
k_values <- 0:9
# Calcolo delle probabilità per ogni valore in k_values
probabilities <- dpois(k_values, lambda = lam_true)
# Stampa delle probabilità
for (i in seq_along(k_values)) {
cat(sprintf("Probabilità di %d eventi: %.4f\n", k_values[i], probabilities[i]))
}
#> Probabilità di 0 eventi: 0.1353
#> Probabilità di 1 eventi: 0.2707
#> Probabilità di 2 eventi: 0.2707
#> Probabilità di 3 eventi: 0.1804
#> Probabilità di 4 eventi: 0.0902
#> Probabilità di 5 eventi: 0.0361
#> Probabilità di 6 eventi: 0.0120
#> Probabilità di 7 eventi: 0.0034
#> Probabilità di 8 eventi: 0.0009
#> Probabilità di 9 eventi: 0.0002Il seguente codice R genera il grafico della funzione di massa di probabilità (PMF) di una distribuzione di Poisson con parametro \(\lambda = 2\):
# Definiamo il parametro lambda
lambd <- 2
# Generiamo i valori sull'asse x (numero di eventi)
x <- 0:9
# Calcoliamo le probabilità corrispondenti
y <- dpois(x, lambda = lambd)
# Creare un dataframe contenente i dati per il grafico
df <- data.frame(
numero_eventi = x,
probabilita = y
)
# Creare il grafico
ggplot(df, aes(x = numero_eventi, y = probabilita)) +
geom_segment(aes(x = numero_eventi, xend = numero_eventi, y = 0, yend = probabilita),
color = "blue", size = 1) + # Segmenti verticali
geom_point(color = "red", size = 3) + # Cerchi sulla cima dei segmenti
labs(title = "Distribuzione di Poisson (λ = 2)",
x = "Numero di eventi",
y = "Probabilità") +
theme(plot.title = element_text(hjust = 0.5)) # Centrare il titolo
52.3 Distribuzione Gamma
La distribuzione Gamma riveste un ruolo centrale nel modello coniugato Gamma-Poisson, poiché funge da distribuzione a priori per il parametro di tasso \(\lambda\) della distribuzione di Poisson. La sua scelta è motivata dalla proprietà di coniugatezza, che consente di ottenere una distribuzione a posteriori appartenente alla stessa famiglia della prior, semplificando i calcoli inferenziali e aggiornando in modo diretto le credenze alla luce dei dati osservati.
La funzione di densità della distribuzione Gamma è definita come:
\[ f(x \mid \alpha, \beta) = \frac{\beta^\alpha}{\Gamma(\alpha)} x^{\alpha - 1} e^{-\beta x}, \]
dove:
- \(\alpha\) (detto parametro di forma) determina la forma della distribuzione: valori più elevati di \(\alpha\) tendono a rendere la distribuzione più simmetrica;
- \(\beta\) (detto parametro di tasso o scala inversa) controlla la scala: valori più elevati di \(\beta\) concentrano maggiormente la massa di probabilità vicino all’origine.
Ad esempio:
- Una distribuzione Gamma con \(\alpha = 2\) e \(\beta = 3\) rappresenta un processo in cui eventi relativamente rari si verificano occasionalmente.
- Una distribuzione Gamma con \(\alpha = 10\) e \(\beta = 1\) descrive un processo più regolare, con eventi che si verificano con maggiore frequenza e prevedibilità.
Per calcolare la densità di probabilità in R, si utilizza:
dgamma(x, shape = alpha, rate = beta)# Definizione dei parametri
alpha <- 2
beta <- 3
# Generazione dei valori di x
x <- seq(0, 3, length.out = 500)
# Calcolo della densità di probabilità
pdf <- dgamma(x, shape = alpha, rate = beta)
# Creare un dataframe contenente i dati per il grafico
df <- data.frame(
x = x,
densita = pdf
)
# Creare il grafico con ggplot2
ggplot(df, aes(x = x, y = densita)) +
geom_line(color = "blue", size = 1) + # Aggiungere la linea per la densità
labs(
title = expression("Distribuzione Gamma (" ~ alpha == 2 ~ "," ~ beta == 3 ~ ")"), # Titolo con espressioni matematiche
x = "x", # Label dell'asse x
y = "Densità di probabilità" # Label dell'asse y
) +
theme(plot.title = element_text(hjust = 0.5)) # Centrare il titolo
Questo grafico mostra la densità di una distribuzione Gamma con \(\alpha = 2\) e \(\beta = 3\), calcolata su un intervallo \([0, 3]\).
52.4 Metodo Basato su Griglia
Vogliamo calcolare la distribuzione a posteriori del parametro \(\lambda\) di un modello di Poisson, usando una distribuzione a priori Gamma. Di seguito mostriamo come si può procedere numericamente, discretizzando lo spazio dei possibili valori di \(\lambda\).
Dati e Assunzioni.
I dati osservati sono:
y <- c(2, 1, 3, 2, 2, 1, 1, 1)e la nostra distribuzione a priori è \(\text{Gamma}(\alpha=9, \beta=2)\).
1. Definizione della Griglia e della Distribuzione A Priori.
- Creiamo una griglia di valori di \(\lambda\), per esempio tra 0.01 e 10, con 1000 punti equidistanti.
- Per ciascun valore di \(\lambda\), calcoliamo la densità della a priori Gamma corrispondente.
È spesso utile visualizzare questa a priori e confrontarla con altre informazioni, ad esempio la media di \(\lambda\) secondo la a priori (\(\alpha/\beta\)) e la media campionaria (\(\bar{y}\)):
# Calcolo di media a priori (media di una Gamma) e media campionaria
media_gamma <- alpha_prior / beta_prior
media_campionaria <- mean(y)
# Creazione di un dataframe per ggplot
df_prior <- data.frame(
lambda = lambda_grid,
densita = prior
)
ggplot(df_prior, aes(x = lambda, y = densita)) +
geom_line() +
geom_vline(xintercept = media_gamma, linetype = "dashed") +
geom_vline(xintercept = media_campionaria, linetype = "dashed") +
labs(
title = "Distribuzione Gamma a Priori",
x = expression(lambda),
y = "Densità"
)
2. Calcolo della Verosimiglianza.
Per una famiglia di distribuzioni di Poisson, la verosimiglianza di un insieme di dati \(y_1, y_2, \dots, y_n\) è il prodotto delle probabilità:
\[ \mathcal{L}(\lambda) \;=\; \prod_{i=1}^n \mathrm{P}(Y=y_i \mid \lambda) \;=\; \prod_{i=1}^n \mathrm{dpois}(y_i, \lambda). \]
In R, possiamo calcolarla per tutti i valori di \(\lambda\) della griglia in modo vettoriale:
L’istruzione precedente fa questo:
-
Passo su ciascun valore di \(\lambda\) nella griglia (
lambda_grid).
- Per ogni \(\lambda\), calcola \(\mathrm{dpois}(y, \lambda)\), che restituisce un vettore di lunghezza \(n\) (il numero di osservazioni), ossia la probabilità di osservare ciascun \(y_i\) secondo una Poisson(\(\lambda\)).
-
Moltiplica insieme tutte queste probabilità con
prod(...), ottenendo un singolo numero che è la verosimiglianza \(\mathcal{L}(\lambda)\) per quell’insieme di dati.
- Ripete il procedimento per ogni \(\lambda\) della griglia, restituendo un vettore di verosimiglianze (una per ciascun punto della griglia).
In altre parole, ogni elemento del vettore likelihood corrisponde alla verosimiglianza calcolata in un punto \(\lambda\). E la verosimiglianza è appunto il prodotto delle probabilità di ciascuno dei valori \(y_i\), assumendo che i dati provengano da una \(\mathrm{Poisson}(\lambda)\) e che le osservazioni siano indipendenti.
3. Calcolo della A Posteriori (non normalizzata).
La distribuzione a posteriori, a meno di una costante di normalizzazione, è data da:
\[ \pi(\lambda \mid y) \; \propto \; \mathcal{L}(\lambda)\,\pi(\lambda). \]
Tradotto in codice:
posterior_unnormalized <- prior * likelihood4. Normalizzazione.
Per ottenere la vera distribuzione di probabilità, dobbiamo normalizzare:
\[ \pi(\lambda \mid y) = \frac{\mathcal{L}(\lambda)\,\pi(\lambda)}{\int \mathcal{L}(\lambda)\,\pi(\lambda)\, d\lambda}. \]
Sul computer, l’integrale si trasforma in una somma pesata dallo “spessore” di ogni intervallo nella griglia:
# Calcolo del fattore di normalizzazione
grid_width <- lambda_grid[2] - lambda_grid[1] # Ampiezza della griglia (costante)
normalization_factor <- sum(posterior_unnormalized) * grid_width
# Distribuzione a posteriori normalizzata
posterior <- posterior_unnormalized / normalization_factor5. Confronto Visivo: A Priori vs A Posteriori.
Infine, possiamo confrontare le due curve (a priori e a posteriori) su un unico grafico:
# Creazione del dataframe
df <- data.frame(
lambda = lambda_grid,
prior = prior,
posterior = posterior
)
# Trasformiamo in formato "lungo" per ggplot
library(reshape2)
df_long <- melt(
df,
id.vars = "lambda",
measure.vars = c("prior", "posterior"),
variable.name = "Distribuzione",
value.name = "Densita"
)
# Grafico
ggplot(df_long, aes(x = lambda, y = Densita, color = Distribuzione)) +
geom_line() +
labs(
title = "Distribuzione A Priori e A Posteriori",
x = expression(lambda),
y = "Densità"
)
Interpretazione.
- La distribuzione a posteriori risulta spostata verso valori più bassi rispetto alla a priori. Questo indica che i dati osservati suggeriscono \(\lambda\) più piccolo di quanto ipotizzato inizialmente dalla a priori.
- La a posteriori appare più concentrata (più piccata) rispetto alla a priori. Ciò riflette la minor incertezza residua dopo aver incorporato l’informazione fornita dai dati.
Questo semplice approccio su griglia mostra chiaramente come la verosimiglianza, “pesata” dalla a priori, generi la nuova distribuzione a posteriori. È un metodo intuitivo che permette di visualizzare passo dopo passo in che modo i dati aggiornano le nostre credenze sul parametro \(\lambda\).
52.5 Modello Coniugato Gamma-Poission
Per calcolare analiticamente la distribuzione a posteriori nel contesto di un modello gamma-Poisson possiamo seguire un processo diretto. Il modello Gamma-Poisson è coniugato, il che significa che la distribuzione a posteriori sarà ancora una distribuzione Gamma.
Seguendo il teorema di Bayes, possiamo scrivere la distribuzione a posteriori come:
\(f(\lambda \mid y) \propto f(y \mid \lambda) \cdot f(\lambda) ,\)
dove \(f(\lambda \mid y)\) è la distribuzione a posteriori, \(f(y \mid \lambda)\) è la verosimiglianza e \(f(\lambda)\) è la distribuzione a priori.
Definiamo la verosimiglianza (distribuzione di Poisson):
\(f(y \mid \lambda) = \prod_{i=1}^n \frac{e^{-\lambda}\lambda^{y_i}}{y_i!} = \frac{e^{-n\lambda}\lambda^{\sum_{i=1}^n y_i}}{\prod_{i=1}^n y_i!}.\)
Definiamo la distribuzione a priori (distribuzione Gamma):
\(f(\lambda) = \frac{b^a}{\Gamma(a)}\lambda^{a-1}e^{-b\lambda}.\)
Ora, moltiplichiamo la verosimiglianza per la distribuzione a priori:
\(f(\lambda|y) \propto \frac{e^{-n\lambda}\lambda^{\sum_{i=1}^n y_i}}{\prod_{i=1}^n y_i!} \cdot \frac{b^a}{\Gamma(a)}\lambda^{a-1}e^{-b\lambda}.\)
Semplifichiamo, eliminando i termini costanti (che non dipendono da \(\lambda\)):
\(f(\lambda \mid y) \propto e^{-n\lambda}\lambda^{\sum_{i=1}^n y_i} \cdot \lambda^{a-1}e^{-b\lambda}.\)
Raggruppiamo i termini:
\(f(\lambda|y) \propto \lambda^{\sum_{i=1}^n y_i} \cdot \lambda^{a-1} \cdot e^{-n\lambda} \cdot e^{-b\lambda}.\)
Semplifichiamo ulteriormente:
\(f(\lambda \mid y) \propto \lambda^{\sum_{i=1}^n y_i + a - 1} \cdot e^{-(n+b)\lambda}.\)
Riconosciamo che questa è la forma di una distribuzione Gamma con nuovi parametri:
\(f(\lambda \mid y) \propto \lambda^{(\sum_{i=1}^n y_i + a) - 1} \cdot e^{-(n+b)\lambda}.\)
Quindi, la distribuzione a posteriori è una Gamma con parametri:
- \(\alpha_{post} = a + \sum_{i=1}^n y_i\),
- \(\beta_{post} = b + n\),
dove:
- \(a\) e \(b\) sono i parametri della distribuzione Gamma a priori,
- \(\sum_{i=1}^n y_i\) è la somma di tutte le osservazioni,
- \(n\) è il numero di osservazioni.
Questa derivazione mostra come la distribuzione a posteriori mantiene la forma di una Gamma, ma con parametri aggiornati che incorporano l’informazione dai dati osservati.
Consideriamo nuovamente l’esempio precedente. Utilizzando i parametri aggiornati, rappresentiamo graficamente la distribuzione a posteriori:
# Parametri aggiornati della distribuzione a posteriori Gamma
alpha_post <- alpha_prior + sum(y)
beta_post <- beta_prior + length(y)
# Griglia dei valori di lambda
lambda_grid <- seq(0.01, 10, length.out = 1000)
# Distribuzione a posteriori analitica
posterior_analytic <- dgamma(lambda_grid, shape = alpha_post, rate = beta_post)
# Creare un dataframe contenente i dati per il grafico
df <- data.frame(
lambda = lambda_grid,
densita = posterior_analytic
)
# Calcolo della media a posteriori
media_posterior <- alpha_post / beta_post
# Creare il grafico con ggplot2
ggplot(df, aes(x = lambda, y = densita)) +
geom_line(color = "blue", size = 1) + # Aggiungere la linea per la densità a posteriori
geom_vline(xintercept = media_posterior, color = "red", linetype = "dashed", size = 1) + # Linea verticale per la media a posteriori
labs(
title = "Distribuzione a Posteriori Analitica Gamma-Poisson",
x = expression(lambda), # Label dell'asse x usando espressioni matematiche
y = "Densità di probabilità" # Label dell'asse y
) +
theme(plot.title = element_text(hjust = 0.5)) + # Centrare il titolo
annotate("text", x = media_posterior + 0.5, y = max(posterior_analytic) * 0.8,
label = "Media a Posteriori", color = "red", size = 4) + # Etichetta per la media a posteriori
scale_x_continuous(expand = expansion(mult = c(0, 0.1))) # Aumentare lo spazio sull'asse x per migliorare la visualizzazione
Il grafico mostra la distribuzione a posteriori analitica del parametro di tasso \(\lambda\) di un modello di Poisson, ottenuta utilizzando una distribuzione a priori Gamma e aggiornando i parametri alla luce dei dati osservati. La distribuzione a posteriori è calcolata come una Gamma aggiornata con i parametri \(\alpha_{\text{post}}\) e \(\beta_{\text{post}}\), e rappresenta la nostra conoscenza aggiornata dopo aver visto i dati. I risultati analitici concordano con quelli ottenuti tramite simulazione.
Procediamo ora con il calcolo della soluzione analitica per la media della distribuzione a posteriori del parametro \(\lambda\):
# Media della distribuzione a posteriori
posterior_mean <- alpha_post / beta_post
cat(sprintf("Media a Posteriori = %.3f\n", posterior_mean))
#> Media a Posteriori = 2.200
# Parametri della distribuzione a posteriori
cat(sprintf("Shape (α) = %.1f\n", alpha_post))
#> Shape (α) = 22.0
cat(sprintf("Rate (β) = %.1f\n", beta_post))
#> Rate (β) = 10.0Possiamo calcolare la probabilità di qualsiasi evento di interesse. Per esempio, ci possiamo chiedere quale sia la probabilità di osservare più di 3 compulsioni per ora:
52.6 Riflessioni Conclusive
Il modello Gamma-Poisson offre un framework robusto per l’inferenza bayesiana su dati di conteggio in psicologia. Partendo da una distribuzione a priori Gamma, che rappresenta la nostra conoscenza iniziale sul tasso medio, siamo in grado di aggiornare questa conoscenza alla luce dei dati osservati, ottenendo una distribuzione a posteriori che riflette in modo più preciso la realtà sottostante. Questo approccio permette di quantificare l’incertezza associata alle nostre stime e di prendere decisioni informate sulla base dei dati disponibili.
52.7 Esercizi
52.8 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] reshape2_1.4.4 mice_3.18.0 thematic_0.1.7 MetBrewer_0.2.0
#> [5] ggokabeito_0.1.0 see_0.11.0 gridExtra_2.3 patchwork_1.3.0
#> [9] bayesplot_1.13.0 psych_2.5.3 scales_1.4.0 markdown_2.0
#> [13] knitr_1.50 lubridate_1.9.4 forcats_1.0.0 stringr_1.5.1
#> [17] dplyr_1.1.4 purrr_1.0.4 readr_2.1.5 tidyr_1.3.1
#> [21] tibble_3.3.0 ggplot2_3.5.2 tidyverse_2.0.0 rio_1.2.3
#> [25] here_1.0.1
#>
#> loaded via a namespace (and not attached):
#> [1] tidyselect_1.2.1 farver_2.1.2 fastmap_1.2.0
#> [4] pacman_0.5.1 digest_0.6.37 rpart_4.1.24
#> [7] timechange_0.3.0 lifecycle_1.0.4 survival_3.8-3
#> [10] magrittr_2.0.3 compiler_4.5.0 rlang_1.1.6
#> [13] tools_4.5.0 labeling_0.4.3 htmlwidgets_1.6.4
#> [16] mnormt_2.1.1 plyr_1.8.9 RColorBrewer_1.1-3
#> [19] withr_3.0.2 nnet_7.3-20 grid_4.5.0
#> [22] jomo_2.7-6 iterators_1.0.14 MASS_7.3-65
#> [25] cli_3.6.5 rmarkdown_2.29 reformulas_0.4.1
#> [28] generics_0.1.4 rstudioapi_0.17.1 tzdb_0.5.0
#> [31] minqa_1.2.8 splines_4.5.0 parallel_4.5.0
#> [34] vctrs_0.6.5 boot_1.3-31 glmnet_4.1-9
#> [37] Matrix_1.7-3 jsonlite_2.0.0 hms_1.1.3
#> [40] mitml_0.4-5 foreach_1.5.2 glue_1.8.0
#> [43] nloptr_2.2.1 pan_1.9 codetools_0.2-20
#> [46] stringi_1.8.7 shape_1.4.6.1 gtable_0.3.6
#> [49] lme4_1.1-37 pillar_1.10.2 htmltools_0.5.8.1
#> [52] R6_2.6.1 Rdpack_2.6.4 rprojroot_2.0.4
#> [55] evaluate_1.0.4 lattice_0.22-7 rbibutils_2.3
#> [58] backports_1.5.0 broom_1.0.8 Rcpp_1.0.14
#> [61] nlme_3.1-168 xfun_0.52 pkgconfig_2.0.3Bibliografia
Scholz, U., Stadler, G., Berli, C., Lüscher, J., & Knoll, N. (2021). How do people experience and respond to social control from their partner? Three daily diary studies. Frontiers in Psychology, 11, 613546.