here::here("code", "_common.R") |>
source()
# Load packages
if (!requireNamespace("pacman")) install.packages("pacman")
pacman::p_load(tidyr)
57 Modello gamma-esponenziale
Prerequisiti
- Leggere il Capitolo 39 e il Capitolo 40 della dispensa.
Concetti e competenze chiave
- Comprendere la distribuzione esponenziale come un modello probabilistico adatto per descrivere i tempi di attesa.
- Sapere applicare il metodo basato su griglia per derivare la distribuzione a posteriori del parametro λ del modello Gamma-Esponenziale.
- Conoscere il modello coniugato Gamma-Esponenziale, dimostrando come la distribuzione a priori Gamma si combini con la verosimiglianza esponenziale per produrre una distribuzione a posteriori Gamma.
- Sapere come calcolare e interpretare le probabilità utilizzando la distribuzione a posteriori.
Preparazione del Notebook
57.1 Introduzione
Nell’inferenza bayesiana, il modello coniugato Gamma-Esponenziale rappresenta un approccio analitico efficace per l’analisi di dati che seguono una distribuzione esponenziale. Questa distribuzione è comunemente utilizzata per modellare i tempi di attesa tra eventi in processi di Poisson, come ad esempio gli intervalli tra arrivi in un sistema a coda o la durata di eventi psicologici.
Consideriamo, ad esempio, un esperimento in psicologia in cui si misurano i tempi di insorgenza di episodi di ansia in seguito a un evento stressante. In questo caso, si può ipotizzare che i tempi di insorgenza siano distribuiti esponenzialmente. Il parametro
Il modello coniugato Gamma-Esponenziale consente di stimare il parametro
57.1.1 Il Modello Matematico
Caso singolo. Supponiamo di osservare un singolo tempo di attesa
dove:
-
rappresenta il tempo di attesa, -
è il parametro della distribuzione, che rappresenta il tasso medio con cui gli eventi si verificano per unità di tempo.
In questa distribuzione,
- Il tempo medio di attesa è il valore medio del tempo che trascorre prima che l’evento si verifichi (ad esempio, quanto tempo ci si aspetta in media prima che arrivi un autobus).
- Il parametro
rappresenta la frequenza con cui l’evento si verifica per unità di tempo, e quindi .
Pertanto, nella funzione esponenziale
Questa funzione descrive la probabilità di osservare un tempo di attesa esattamente uguale a
Per fare un esempio, supponiamo che il tempo medio di attesa sia 2 ore. In questo caso, il parametro
# Parametro lambda (l'inverso del tempo medio di attesa)
lambda_value <- 1 / 2
Disegniamo la funzione di densità esponenziale per tempi di attesa compresi tra 0 e 10 ore.
# Creazione dei valori di y1 su cui valutare la funzione
y1_values <- seq(0, 10, length.out = 500) # Intervallo da 0 a 10 ore
# Definizione della funzione di densità esponenziale
exponential_density <- function(y1, lambda_value) {
lambda_value * exp(-lambda_value * y1)
}
# Calcolo dei valori di f(y1 | lambda)
f_values <- exponential_density(y1_values, lambda_value)
# Creazione del grafico
plot(
y1_values, f_values,
type = "l", lwd = 2,
main = expression(paste(
"Funzione di Densità Esponenziale con ",
lambda, " = 1/2"
)),
xlab = "Tempo di attesa (ore)", ylab = "Densità"
)
legend("topright",
legend = expression(f(y[1] ~ "|" ~ lambda) == lambda ~ e^(-lambda * y[1])),
lty = 1, bty = "n"
)
Questo codice crea un grafico della funzione di densità esponenziale per
Il caso di
Poiché le osservazioni sono indipendenti, la probabilità congiunta di osservare tutti i tempi di attesa nel campione è il prodotto delle probabilità individuali. Di conseguenza, la funzione di verosimiglianza per l’intero campione è data da:
Sostituendo la funzione di densità della distribuzione esponenziale per ciascuna osservazione, otteniamo:
Raccogliendo i termini comuni, possiamo riscrivere la funzione di verosimiglianza come:
Utilizzando la notazione di sommatoria, possiamo esprimere la funzione di verosimiglianza in modo compatto come:
La funzione di verosimiglianza
Spesso è più conveniente lavorare con il logaritmo della funzione di verosimiglianza, poiché il logaritmo trasforma i prodotti in somme, semplificando i calcoli. Il logaritmo della funzione di verosimiglianza è:
Questa forma semplificata della log-verosimiglianza è utile per stimare il parametro
57.2 Aggiornare le Nostre Credenze con l’Inferenza Bayesiana
Nell’approccio bayesiano, non consideriamo solo i dati, ma anche le nostre conoscenze a priori sul parametro
Per fare un esempio concreto, simuleremo un campione di dati. Immaginiamo di raccogliere dati che rappresentano il tempo, misurato in ore, che intercorre tra episodi di ansia in individui con disturbi d’ansia. La distribuzione esponenziale può essere utilizzata per modellare questo tempo di attesa tra un episodio di ansia e il successivo. Ad esempio, possiamo ipotizzare che, una volta concluso un episodio, l’insorgenza del prossimo segua un processo stocastico con un tasso costante, indipendentemente dal tempo trascorso dall’episodio precedente. In questo contesto, il parametro
Se
# Impostiamo il seed per rendere i risultati riproducibili
set.seed(42)
# Parametro medio (tempo medio di attesa in ore)
mean_time <- 3.0
lambda_param <- 1 / mean_time # Parametro lambda per la distribuzione esponenziale
# Numero di osservazioni nel campione
n <- 15
# Generazione del campione casuale
y <- round(rexp(n, rate = lambda_param))
y
#> [1] 1 2 1 0 1 4 1 1 4 2 4 7 0 0 4
Immaginiamo che questi dati rappresentino il tempo di attesa in ore tra episodi di ansia in 15 individui con disturbi d’ansia. Il tempo di attesa medio è:
mean_time_observed <- mean(y)
mean_time_observed
#> [1] 2.133
Il tasso di occorrenza,
lambda_estimated <- 1 / mean_time_observed
lambda_estimated
#> [1] 0.4688
Ora possiamo procedere a trovare la distribuzione a posteriori per il tasso di occorrenza
57.2.1 Passi per definire un prior debolmente informativo
Il primo passo consiste nel definire una distribuzione a priori per
La distribuzione a priori coniugata per la distribuzione esponenziale è la distribuzione Gamma. Un prior debolmente informativo per
Supponiamo di voler impostare una distribuzione a priori per il tasso di occorrenza
I parametri della distribuzione Gamma,
Verifichiamo che questi parametri corrispondano alla media e alla varianza desiderate.
# Parametri della distribuzione Gamma
alpha_prior <- 0.45 # Forma
beta_prior <- 1.5 # Tasso (inverso della scala)
# Calcolo della media e della varianza della distribuzione Gamma per λ
mean_lambda <- alpha_prior / beta_prior
variance_lambda <- alpha_prior / (beta_prior^2)
# Calcolo della media del tempo di attesa (E[Y] = 1 / E[λ])
mean_waiting_time <- 1 / mean_lambda
# Stampiamo la media e la varianza di λ e del tempo di attesa
cat("Media di λ:", mean_lambda, "\n")
#> Media di λ: 0.3
cat("Varianza di λ:", variance_lambda, "\n")
#> Varianza di λ: 0.2
cat("Media del tempo di attesa:", mean_waiting_time, "ore\n")
#> Media del tempo di attesa: 3.333 ore
La media del tempo di attesa è 3.33 ore, come desiderato. La varianza del tempo di attesa richiede un calcolo più complesso e non è semplicemente l’inverso della varianza di
57.2.2 Visualizzazione della Distribuzione a Priori
Disegniamo la funzione di densità della distribuzione Gamma per
# Creazione dei valori x per il grafico
x <- seq(0, 2, length.out = 500)
# Funzione di densità della distribuzione Gamma con i parametri dati
gamma_pdf <- dgamma(x, shape = alpha_prior, rate = beta_prior)
# Creazione del grafico
plot(x, gamma_pdf,
type = "l", lwd = 2,
main = bquote("Funzione di Densità Gamma con " ~ alpha == .(alpha_prior) ~ "," ~ beta == .(beta_prior)),
xlab = expression(lambda),
ylab = "Densità di probabilità"
)
legend(
"topright",
legend = "Distribuzione a Priori",
lty = 1,
lwd = 2,
bty = "n"
)
57.3 Metodo Basato su Griglia
Poniamoci l’obiettivo di utilizzare il metodo basato su griglia per derivare la distribuzione a posteriori del parametro
# Evitiamo zero per evitare divisione per zero
lambda_grid <- seq(0.01, 2, length.out = 1000)
57.3.1 Calcolo della Distribuzione a Priori
# Calcolo della distribuzione a priori
prior <- dgamma(lambda_grid, shape = alpha_prior, rate = beta_prior)
57.3.2 Calcolo della Verosimiglianza
Immaginiamo di avere i seguenti dati, che rappresentano i tempi di attesa (in ore) tra episodi di ansia per 15 individui:
# Impostiamo il seed per rendere i risultati riproducibili
set.seed(42)
# Parametro medio (tempo medio di attesa in ore)
mean_time <- 3.0
lambda_param <- 1 / mean_time # Parametro λ per la distribuzione esponenziale
# Numero di osservazioni nel campione
n <- 15
# Generazione del campione casuale
y <- round(rexp(n, rate = lambda_param))
y
#> [1] 1 2 1 0 1 4 1 1 4 2 4 7 0 0 4
Calcoliamo la log-verosimiglianza per ciascun valore di
57.3.3 Calcolo della Distribuzione a Posteriori
Calcoliamo la distribuzione a posteriori non normalizzata in logaritmo per evitare problemi di underflow numerico:
# Calcolo del log-posteriori non normalizzato
log_posterior_unnormalized <- log_likelihood + log(prior)
# Normalizzazione per stabilità numerica
log_posterior_unnormalized <- log_posterior_unnormalized - max(log_posterior_unnormalized)
# Convertiamo in scala normale
posterior_unnormalized <- exp(log_posterior_unnormalized)
Normalizziamo la distribuzione a posteriori:
57.3.4 Visualizzazione dei Risultati
# Grafico della distribuzione a posteriori e a priori
plot(lambda_grid, posterior,
type = "l",
lwd = 2,
xlab = expression(lambda),
ylab = "Densità di probabilità",
main = "Distribuzione a Posteriori di λ"
)
lines(lambda_grid, prior, lwd = 2, lty = 2)
legend("topright",
legend = c("Posteriori", "Priori"),
lty = c(1, 2),
lwd = 2,
bty = "n"
)
Dal grafico, possiamo osservare come la distribuzione a posteriori di
57.3.5 Calcolo della Media e Varianza Posteriori
Possiamo calcolare la media e la varianza a posteriori di
# Calcolo della media a posteriori di λ
posterior_mean_lambda <- sum(lambda_grid * posterior * delta_lambda)
# Calcolo della varianza a posteriori di λ
posterior_variance_lambda <-
sum((lambda_grid - posterior_mean_lambda)^2 * posterior * delta_lambda)
# Stampiamo i risultati
cat("Media a posteriori di λ:", posterior_mean_lambda, "\n")
#> Media a posteriori di λ: 0.4612
cat("Varianza a posteriori di λ:", posterior_variance_lambda, "\n")
#> Varianza a posteriori di λ: 0.01377
57.3.6 Stima del Tempo di Attesa Medio a Posteriori
Utilizzando la media a posteriori di
# Stima del tempo di attesa medio a posteriori
posterior_mean_waiting_time <- 1 / posterior_mean_lambda
cat("Tempo di attesa medio a posteriori:", posterior_mean_waiting_time, "ore\n")
#> Tempo di attesa medio a posteriori: 2.168 ore
In conclusione, attraverso il metodo basato su griglia, abbiamo derivato la distribuzione a posteriori del tasso di occorrenza
57.4 Modello Coniugato Gamma-Esponenziale
Quando utilizziamo una distribuzione
In altre parole, se il parametro
Questo aggiornamento dei parametri è una conseguenza della proprietà coniugata della distribuzione Gamma rispetto alla distribuzione esponenziale.
57.4.1 Dimostrazione del Modello Coniugato Gamma-Esponenziale
Per dimostrare questo risultato, partiamo dal teorema di Bayes:
dove
La funzione di verosimiglianza per un campione di
Supponiamo che il parametro
Moltiplicando la verosimiglianza per la distribuzione a priori, otteniamo la distribuzione a posteriori:
Semplificando, si ottiene:
Questa espressione corrisponde alla forma di una distribuzione Gamma con parametri aggiornati:
- parametro della forma (alpha):
; - parametro della scala (beta):
.
Quindi, la distribuzione a posteriori di
Questa derivazione mostra come l’informazione contenuta nei dati osservati venga incorporata nei parametri della distribuzione a posteriori, mantenendo la forma della distribuzione a priori grazie alla proprietà coniugata.
57.4.2 Calcolo dei Parametri Aggiornati
Per il caso dell’esempio in discussione:
- Il numero di osservazioni nel campione
è 15; - La somma delle osservazioni nel campione è:
La distribuzione a posteriori per
Usando i valori
# Dati iniziali
alpha_prior <- 0.009
beta_prior <- 0.03
n <- 15
y <- c(1, 9, 4, 3, 1, 1, 0, 6, 3, 4, 0, 11, 5, 1, 1)
sum_y <- sum(y)
# Parametri aggiornati
alpha_post <- alpha_prior + n
beta_post <- beta_prior + sum_y
# Stampa dei parametri aggiornati
cat(sprintf("alpha_post = %.3f; beta_post = %.3f\n", alpha_post, beta_post))
#> alpha_post = 15.009; beta_post = 50.030
Generiamo il grafico della distribuzione a posteriori.
# Griglia di valori di lambda per il grafico
lambda_grid <- seq(0, 2, length.out = 1000)
# Calcolo della densità della distribuzione Gamma a posteriori
posterior_pdf <- dgamma(lambda_grid, shape = alpha_post, rate = beta_post)
# Creazione del grafico
plot(
lambda_grid, posterior_pdf,
type = "l", lwd = 2,
main = bquote(
"Distribuzione Gamma a Posteriori con " ~ alpha[post] ==
.(alpha_post) ~ "," ~ beta[post] == .(beta_post)
),
xlab = expression(lambda), ylab = "Densità di probabilità"
)
legend(
"topright",
legend = bquote(
Gamma(alpha[post] == .(alpha_post) ~ "," ~ beta[post] == .(beta_post))
),
lty = 1,
lwd = 2,
bty = "n"
)
57.4.3 Calcolo della Media e della Varianza a Posteriori
La media e la varianza della distribuzione a posteriori sono calcolate come:
In R:
# Calcolo della media e della varianza a posteriori
posterior_mean <- alpha_post / beta_post
posterior_variance <- alpha_post / (beta_post^2)
# Stampa dei risultati
cat(sprintf("Media a posteriori di λ: %.3f\n", posterior_mean))
#> Media a posteriori di λ: 0.300
cat(sprintf("Varianza a posteriori di λ: %.6f\n", posterior_variance))
#> Varianza a posteriori di λ: 0.005996
La distribuzione a posteriori di
-
Media a Posteriori: La media a posteriori di
è circa 0.3, indicando un tasso di occorrenza medio di 0.3 episodi per ora. -
Varianza a Posteriori: La varianza di 0.006 riflette una moderata incertezza nella stima di
.
Questi risultati evidenziano come la distribuzione Gamma aggiornata combini informazioni a priori e dati osservati per ottenere stime robuste del tasso di occorrenza
57.4.4 Trasformazione della media e della varianza
Per interpretare questi risultati in termini di tempi di attesa, dobbiamo trasformare i valori sulla scala inversa, cioè passare dal tasso
La media del tempo di attesa
Per la varianza del tempo di attesa, dobbiamo applicare una trasformazione più complessa. La varianza del tempo di attesa
Sostituendo i valori calcolati per la varianza e la media di
In conclusione,
- la media del tempo di attesa a posteriori è di circa 3.33 ore;
- la varianza del tempo di attesa a posteriori è di circa 24.59 ore, che indica una certa dispersione dei tempi di attesa, con alcuni episodi che possono verificarsi molto più rapidamente o più lentamente rispetto alla media. Si noti che la distribuzione a posteriori è più concentrata rispetto al prior, poiché incorpora l’informazione aggiuntiva proveniente dai dati osservati.
Questa trasformazione ci permette di interpretare i risultati sulla scala dei tempi di attesa, che è più intuitiva per descrivere fenomeni psicologici come l’insorgenza di episodi di ansia.
57.5 Applicazioni
Una volta ottenuta la distribuzione a posteriori per λ, possiamo utilizzarla per rispondere a domande probabilistiche relative ai tempi di attesa tra episodi di ansia. Questo approccio ci consente di calcolare probabilità aggiornate alla luce dei dati osservati, rispecchiando meglio l’incertezza e le informazioni disponibili. Ad esempio, possiamo stimare la probabilità di osservare tempi di attesa compresi tra 2 e 5 ore.
Per risolvere questo problema possiamo usare il metodo Monte Carlo:
- Generiamo un gran numero di campioni dalla distribuzione Gamma a posteriori di
, usando i parametri e . - Convertiamo ciascun
campionato in un tempo di attesa . - Calcoliamo le probabilità richieste (ad esempio,
e ) semplicemente contando le proporzioni di campioni che soddisfano tali condizioni.
57.5.1 Simulazione con Monte Carlo
Utilizziamo i parametri posteriori della distribuzione Gamma per simulare campioni di
# Parametri posteriori per la distribuzione Gamma
alpha_post <- 15.009
beta_post <- 50.03
# Numero di campioni da generare
n_samples <- 100000
# Simulazione di campioni dalla distribuzione Gamma per λ
lambda_samples <- rgamma(n_samples, shape = alpha_post, rate = beta_post)
# Conversione dei campioni di λ in tempi di attesa T = 1/λ
waiting_time_samples <- 1 / lambda_samples
# Calcolo della probabilità che il tempo di attesa sia compreso tra 2 e 5 ore
prob_between_2_and_5_mc <-
mean(waiting_time_samples >= 2 & waiting_time_samples <= 5)
# Stampa del risultato
cat(sprintf("Probabilità che il tempo di attesa sia tra 2 e 5 ore: %.4f\n", prob_between_2_and_5_mc))
#> Probabilità che il tempo di attesa sia tra 2 e 5 ore: 0.9027
La probabilità che il tempo di attesa tra due episodi di ansia sia compreso tra 2 e 5 ore è di circa 0.9029, ovvero il 90.3%. Questa stima si basa su campioni simulati dalla distribuzione a posteriori di
57.6 Riflessioni Conclusive
Il modello esponenziale si rivela uno strumento utile e versatile nella ricerca psicologica, in particolare per la modellazione di fenomeni caratterizzati da tempi di attesa o durate di eventi. È applicabile in vari contesti, tra cui l’analisi dei tempi di reazione, lo studio degli intervalli tra episodi di ansia o depressione e, più in generale, in tutti quei processi psicologici che possono essere descritti come il tempo trascorso fino al verificarsi di un evento.
Un aspetto particolarmente vantaggioso dell’uso di modelli basati sulla distribuzione esponenziale nell’inferenza bayesiana è la possibilità di utilizzare le famiglie coniugate. Nella famiglia coniugata Gamma-Esponenziale, la distribuzione a priori Gamma per il tasso
La combinazione tra la distribuzione esponenziale e il prior Gamma non solo semplifica l’inferenza, ma fornisce anche una struttura interpretativa chiara. Ad esempio, il parametro
Inoltre, grazie alla flessibilità del metodo Monte Carlo, è possibile simulare campioni dalla distribuzione a posteriori di
In conclusione, il modello esponenziale, integrato in un framework bayesiano con famiglie coniugate, rappresenta un utile strumento per la ricerca psicologica. Offre un modo rigoroso per modellare e analizzare dati su tempi di attesa e durate, fornendo al contempo una base solida per l’inferenza e la previsione dei processi psicologici.
Informazioni sull’Ambiente di Sviluppo
sessionInfo()
#> R version 4.4.2 (2024-10-31)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS Sequoia 15.3.2
#>
#> Matrix products: default
#> BLAS: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRblas.0.dylib
#> LAPACK: /Library/Frameworks/R.framework/Versions/4.4-arm64/Resources/lib/libRlapack.dylib; LAPACK version 3.12.0
#>
#> 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] thematic_0.1.6 MetBrewer_0.2.0 ggokabeito_0.1.0 see_0.11.0
#> [5] gridExtra_2.3 patchwork_1.3.0 bayesplot_1.11.1 psych_2.5.3
#> [9] scales_1.3.0 markdown_2.0 knitr_1.50 lubridate_1.9.4
#> [13] forcats_1.0.0 stringr_1.5.1 dplyr_1.1.4 purrr_1.0.4
#> [17] readr_2.1.5 tidyr_1.3.1 tibble_3.2.1 ggplot2_3.5.1
#> [21] tidyverse_2.0.0 rio_1.2.3 here_1.0.1
#>
#> loaded via a namespace (and not attached):
#> [1] generics_0.1.3 stringi_1.8.4 lattice_0.22-6 hms_1.1.3
#> [5] digest_0.6.37 magrittr_2.0.3 evaluate_1.0.3 grid_4.4.2
#> [9] timechange_0.3.0 fastmap_1.2.0 rprojroot_2.0.4 jsonlite_1.9.1
#> [13] mnormt_2.1.1 cli_3.6.4 rlang_1.1.5 munsell_0.5.1
#> [17] withr_3.0.2 tools_4.4.2 parallel_4.4.2 tzdb_0.5.0
#> [21] colorspace_2.1-1 pacman_0.5.1 vctrs_0.6.5 R6_2.6.1
#> [25] lifecycle_1.0.4 htmlwidgets_1.6.4 pkgconfig_2.0.3 pillar_1.10.1
#> [29] gtable_0.3.6 glue_1.8.0 xfun_0.51 tidyselect_1.2.1
#> [33] rstudioapi_0.17.1 farver_2.1.2 htmltools_0.5.8.1 nlme_3.1-167
#> [37] rmarkdown_2.29 compiler_4.4.2