54  Modello coniugato Gamma-Poisson

In questo capitolo imparerai a
  • Comprendere la distribuzione di Poisson come un modello probabilistico adatto per descrivere eventi rari in un intervallo di tempo o spazio fisso.
  • Sapere applicare il metodo basato su griglia per derivare la distribuzione a posteriori del parametro \(\lambda\) della distribuzione di Poisson.
  • Conoscere il modello coniugato Gamma-Poisson, dimostrando come la distribuzione a priori Gamma si combini con la verosimiglianza di Poisson per produrre una distribuzione a posteriori Gamma.
  • Sapere come calcolare e interpretare le probabilità utilizzando la distribuzione a posteriori.
Prerequisiti
Preparazione del Notebook 
here::here("code", "_common.R") |> 
  source()

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

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

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

54.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.0002

Il 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

54.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à.

In R, la parametrizzazione della distribuzione Gamma utilizza direttamente il parametro \(\beta\) come tasso (\(rate\)), che è l’inverso del parametro di scala (\(scale = 1 / \beta\)).

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

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

  1. Creiamo una griglia di valori di \(\lambda\), per esempio tra 0.01 e 10, con 1000 punti equidistanti.
  2. Per ciascun valore di \(\lambda\), calcoliamo la densità della a priori Gamma corrispondente.
# Parametri della distribuzione a priori Gamma
alpha_prior <- 9
beta_prior  <- 2

# Griglia di valori di lambda
lambda_grid <- seq(0.01, 10, length.out = 1000)

# Densità della distribuzione a priori
prior <- dgamma(lambda_grid, shape = alpha_prior, rate = beta_prior)

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

# Calcolo vettoriale della verosimiglianza
likelihood <- sapply(lambda_grid, function(l) prod(dpois(y, l)))

L’istruzione precedente fa questo:

  1. Passo su ciascun valore di \(\lambda\) nella griglia (lambda_grid).
  2. 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\)).
  3. Moltiplica insieme tutte queste probabilità con prod(...), ottenendo un singolo numero che è la verosimiglianza \(\mathcal{L}(\lambda)\) per quell’insieme di dati.
  4. 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 * likelihood

4. 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_factor

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

  1. 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.
  2. 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\).

54.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.0

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

# Probabilità di osservare più di 3 eventi
prob_y_greater_than_3 <- 1 - pgamma(3, shape = alpha_post, rate = beta_post)
cat(sprintf("Probabilità di osservare più di 3 eventi = %.3f\n", prob_y_greater_than_3))
#> Probabilità di osservare più di 3 eventi = 0.054

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

54.7 Esercizi

Consideriamo uno studio longitudinale su coppie, dove i partecipanti registrano quotidianamente la frequenza con cui nascondono il loro comportamento di fumo al partner. Basandoci sui dati di Scholz et al. (2021), assumiamo che il tasso medio di nascondere il fumo sia di 1.52 volte al giorno. Supponiamo di avere i seguenti dati giornalieri per un partecipante:

  • Giorno 1: 2 volte.
  • Giorno 2: 0 volte.
  • Giorno 3: 1 volta.
  • Giorno 4: 3 volte.

Utilizzare un modello Gamma-Poisson per stimare la distribuzione a posteriori del tasso individuale di nascondere il fumo per un partecipante specifico, dato il suo insieme di osservazioni giornaliere.

Vogliamo stimare quante volte al giorno una persona tende a nascondere il proprio fumo al partner. Abbiamo quattro giorni di dati su questa persona e, grazie a un precedente studio, sappiamo che in media le persone lo fanno circa 1.52 volte al giorno.

Ma i dati di una sola persona sono pochi, quindi useremo un approccio bayesiano per combinare:

  1. Le informazioni preesistenti (dal precedente studio),
  2. Le nuove osservazioni (quanti eventi sono stati registrati nei quattro giorni).

Con un modello Gamma-Poisson, possiamo aggiornare la nostra stima e ottenere una distribuzione a posteriori, che ci dirà quali sono i valori più probabili per il tasso di nascondimento di questa persona.

# Dati osservati: numero di volte che la persona ha nascosto il fumo ogni giorno
osservazioni <- c(2, 0, 1, 3)

# Numero totale di giorni osservati
n_giorni <- length(osservazioni)

# Numero totale di eventi (quante volte ha nascosto il fumo in totale)
eventi_totali <- sum(osservazioni)

# Informazione a priori dallo studio precedente
media_priori <- 1.52
forma_priori <- 3   # Parametro di forma scelto per un'incertezza moderata
tasso_priori <- forma_priori / media_priori  # Parametro di scala

# Aggiornamento bayesiano: parametri della distribuzione a posteriori
forma_post <- forma_priori + eventi_totali
tasso_post <- tasso_priori + n_giorni

# Creazione della griglia di valori possibili per il tasso giornaliero (lambda)
lambda_valori <- seq(0, 5, length.out = 100)

# Calcolo delle densità per le distribuzioni a priori e a posteriori
densita_priori <- dgamma(lambda_valori, shape = forma_priori, rate = tasso_priori)
densita_post <- dgamma(lambda_valori, shape = forma_post, rate = tasso_post)

# Creazione di un dataframe per il grafico
dati_plot <- data.frame(
  lambda = rep(lambda_valori, 2),
  densita = c(densita_priori, densita_post),
  distribuzione = rep(c("Priori", "Posteriori"), each = length(lambda_valori))
)

# Creazione del grafico per confrontare la distribuzione a priori e quella aggiornata (posteriori)
ggplot(dati_plot, aes(x = lambda, y = densita, color = distribuzione)) +
  geom_line(size = 1.2) +
  labs(
    title = "Stima del tasso di nascondimento del fumo",
    x = "Tasso giornaliero (λ)",
    y = "Densità"
  ) +
  theme_minimal() +
  theme(legend.position = "bottom")

Cosa fa questo codice

  1. Carica i dati: abbiamo registrato per 4 giorni quante volte questa persona ha nascosto il fumo.
  2. Imposta una conoscenza iniziale (priori): basata sullo studio precedente.
  3. Applica il teorema di Bayes per aggiornare la stima con i nuovi dati.
  4. Genera il grafico: confronta la distribuzione prima e dopo l’aggiornamento con i dati osservati.

Interpretazione del risultato

  • La curva rossa (priori) rappresenta la nostra stima iniziale, basata sullo studio precedente.
  • La curva blu (posteriori) è la nostra nuova stima dopo aver considerato i dati della persona.
  • Se i dati raccolti sono molto diversi dal valore medio dello studio, la curva blu si sposterà rispetto alla rossa.

Questa analisi ci permette di stimare il comportamento specifico di una persona integrando dati generali e osservazioni individuali, in modo più informativo rispetto a una semplice media.

54.8 Informazioni sull’Ambiente di Sviluppo

sessionInfo()
#> R version 4.4.2 (2024-10-31)
#> Platform: aarch64-apple-darwin20
#> Running under: macOS Sequoia 15.3.1
#> 
#> 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] reshape2_1.4.4   mice_3.17.0      thematic_0.1.6   MetBrewer_0.2.0 
#>  [5] ggokabeito_0.1.0 see_0.10.0       gridExtra_2.3    patchwork_1.3.0 
#>  [9] bayesplot_1.11.1 psych_2.4.12     scales_1.3.0     markdown_1.13   
#> [13] knitr_1.49       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.2.1     ggplot2_3.5.1    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     pacman_0.5.1     
#>  [5] digest_0.6.37     rpart_4.1.24      timechange_0.3.0  lifecycle_1.0.4  
#>  [9] survival_3.8-3    magrittr_2.0.3    compiler_4.4.2    rlang_1.1.5      
#> [13] tools_4.4.2       labeling_0.4.3    htmlwidgets_1.6.4 mnormt_2.1.1     
#> [17] plyr_1.8.9        withr_3.0.2       nnet_7.3-20       grid_4.4.2       
#> [21] jomo_2.7-6        colorspace_2.1-1  iterators_1.0.14  MASS_7.3-65      
#> [25] cli_3.6.4         rmarkdown_2.29    reformulas_0.4.0  generics_0.1.3   
#> [29] rstudioapi_0.17.1 tzdb_0.4.0        minqa_1.2.8       splines_4.4.2    
#> [33] parallel_4.4.2    vctrs_0.6.5       boot_1.3-31       glmnet_4.1-8     
#> [37] Matrix_1.7-2      jsonlite_1.9.1    hms_1.1.3         mitml_0.4-5      
#> [41] foreach_1.5.2     glue_1.8.0        nloptr_2.1.1      pan_1.9          
#> [45] codetools_0.2-20  stringi_1.8.4     shape_1.4.6.1     gtable_0.3.6     
#> [49] lme4_1.1-36       munsell_0.5.1     pillar_1.10.1     htmltools_0.5.8.1
#> [53] R6_2.6.1          Rdpack_2.6.2      rprojroot_2.0.4   evaluate_1.0.3   
#> [57] lattice_0.22-6    rbibutils_2.3     backports_1.5.0   broom_1.0.7      
#> [61] Rcpp_1.0.14       nlme_3.1-167      xfun_0.51         pkgconfig_2.0.3

Bibliografia

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.