41  La verosimiglianza gaussiana

In questo capitolo, imparerai a:
  • comprendere il concetto di verosimiglianza nel contesto del modello gaussiano.
Prerequisiti
Preparazione del Notebook 
here::here("code", "_common.R") |> 
  source()

Introduzione

La distribuzione gaussiana (o distribuzione normale) è una delle distribuzioni più utilizzate in statistica perché descrive molti fenomeni naturali e psicologici. In questo capitolo esploreremo come si calcola la verosimiglianza, ovvero la plausibilità dei parametri, nel caso della distribuzione normale.

41.1 Modello Gaussiano e Verosimiglianza

41.1.1 Caso di una Singola Osservazione

Immaginiamo di misurare il Quoziente Intellettivo (QI) di una persona e ottenere un valore specifico, ad esempio 114. Assumiamo che il QI segua una distribuzione normale con media \(\mu\) sconosciuta e deviazione standard \(\sigma\) nota (ad esempio \(\sigma = 15\)).

La funzione di densità di probabilità per una distribuzione normale è:

\[ f(y \mid \mu, \sigma) = \frac{1}{\sigma\sqrt{2\pi}} \exp\left(-\frac{(y-\mu)^2}{2\sigma^2}\right) , \]

dove:

  • \(y\) è il valore osservato,
  • \(\mu\) è la media (il parametro che vogliamo stimare),
  • \(\sigma\) è la deviazione standard (conosciuta).

La verosimiglianza misura quanto diversi valori di \(\mu\) sono plausibili, dato il valore osservato (114).

Esempio pratico in R:

# Dati iniziali
y_obs <- 114
sigma <- 15
mu_values <- seq(70, 160, length.out = 1000)

# Calcolo della verosimiglianza
likelihood <- dnorm(y_obs, mean = mu_values, sd = sigma)

# Grafico della verosimiglianza
ggplot(data.frame(mu = mu_values, likelihood = likelihood), aes(x = mu, y = likelihood)) +
  geom_line(color = "blue") +
  labs(
    title = "Verosimiglianza per un singolo valore di QI (114)",
    x = "Media (μ)",
    y = "Verosimiglianza"
  )

Qual è il valore migliore per \(\mu\)?

Il valore migliore di \(\mu\) sarà quello che rende massima la verosimiglianza. In questo semplice caso, è esattamente il valore osservato (114):

mu_ottimale <- mu_values[which.max(likelihood)]
cat("Il valore ottimale di μ è:", mu_ottimale)
#> Il valore ottimale di μ è: 114

41.1.2 Log-Verosimiglianza

Spesso, per semplicità di calcolo, si usa la log-verosimiglianza, che trasforma i prodotti in somme, rendendo i calcoli più semplici e stabili:

\[ \log L(\mu \mid y, \sigma) = -\frac{1}{2}\log(2\pi) - \log(\sigma) - \frac{(y-\mu)^2}{2\sigma^2}. \]

Calcolo pratico con R:

negative_log_likelihood <- function(mu, y, sigma) {
  0.5 * log(2 * pi) + log(sigma) + ((y - mu)^2) / (2 * sigma^2)
}

result <- optim(
  par = 100, # Valore iniziale
  fn = negative_log_likelihood,
  y = y_obs,
  sigma = sigma,
  method = "L-BFGS-B",
  lower = 70,
  upper = 160
)

mu_max_loglik <- result$par
cat("Il valore ottimale di μ dalla log-verosimiglianza è:", mu_max_loglik)
#> Il valore ottimale di μ dalla log-verosimiglianza è: 114

In questo caso, otteniamo nuovamente \(\mu = 114\).

41.1.3 Campione di Osservazioni Indipendenti

Consideriamo ora un insieme più grande di osservazioni indipendenti, ciascuna co una distribuzione normale con la stessa media \(\mu\) e la stessa deviazione standard \(\sigma\). Supponiamo di aver misurato i punteggi di depressione (ad esempio, scala BDI-II) per 30 persone:

y <- c(
  26, 35, 30, 25, 44, 30, 33, 43, 22, 43, 24, 19, 39, 31, 25, 28, 35, 30, 26, 
  31, 41, 36, 26, 35, 33, 28, 27, 34, 27, 22
)

Assumiamo inoltre che \(\sigma\) sia nota (ad esempio \(\sigma = 6.50\)). In questo caso, la log-verosimiglianza si calcola così:

log_likelihood <- function(mu, y, sigma) {
  n <- length(y)
  term1 <- -n * log(sigma) - n * log(sqrt(2 * pi))
  term2 <- -sum((y - mu)^2) / (2 * sigma^2)
  return(term1 + term2)
}

sigma <- 6.50
mu_range <- seq(mean(y) - 2 * sigma, mean(y) + 2 * sigma, length.out = 1000)

# Calcolo della log-verosimiglianza con un ciclo for
log_lik_values <- numeric(length(mu_range))

for (i in seq_along(mu_range)) {
  log_lik_values[i] <- log_likelihood(mu_range[i], y, sigma)
}

Il ciclo for() scorre ogni valore di mu_range, calcolando la log-verosimiglianza per ciascun valore e salvando i risultati nel vettore log_lik_values.

Visualizziamo la log-verosimiglianza:

ggplot(
  data.frame(mu = mu_range, log_likelihood = log_lik_values),
  aes(x = mu, y = log_likelihood)
) +
  geom_line(color = "blue") +
  labs(
    title = "Log-verosimiglianza per dati BDI-II",
    x = "Media (μ)",
    y = "Log-verosimiglianza"
  ) +
  geom_vline(xintercept = mean(y), linetype = "dashed", color = "red")

Ottimizzazione numerica per \(\mu\):

negative_log_likelihood <- function(mu, y, sigma) {
  -log_likelihood(mu, y, sigma)
}

result <- optim(
  par = mean(y),
  fn = negative_log_likelihood,
  y = y,
  sigma = sigma,
  method = "L-BFGS-B",
  lower = min(mu_range),
  upper = max(mu_range)
)

mu_optimal <- result$par
cat("Il valore ottimale di μ è:", mu_optimal)
#> Il valore ottimale di μ è: 30.93

41.2 Riflessioni Conclusive

  • Nel caso di una distribuzione normale con \(\sigma\) nota, la miglior stima di \(\mu\) (stima di massima verosimiglianza) è sempre la media delle osservazioni.
  • Visualizzare la funzione di verosimiglianza aiuta a capire la plausibilità relativa dei diversi valori di \(\mu\).
  • La log-verosimiglianza semplifica e rende più stabili i calcoli numerici.

Esercizi

Supponi di aver misurato il livello di ansia (ad esempio usando una scala standardizzata) in un campione di 15 persone con i seguenti punteggi:

ansia <- c(23, 27, 30, 29, 25, 28, 26, 24, 31, 29, 27, 26, 28, 30, 25)

Assumendo che la deviazione standard sia nota e pari a 3.5, svolgi le seguenti attività in R:

  1. Calcola la funzione di verosimiglianza gaussiana per diversi valori di \(\mu\) nell’intervallo da 20 a 35.
  2. Trova numericamente il valore di \(\mu\) che rende massima la verosimiglianza (stima di massima verosimiglianza, MLE).
  3. Disegna un grafico della funzione di verosimiglianza per visualizzare il risultato.

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] thematic_0.1.6   MetBrewer_0.2.0  ggokabeito_0.1.0 see_0.10.0      
#>  [5] gridExtra_2.3    patchwork_1.3.0  bayesplot_1.11.1 psych_2.4.12    
#>  [9] scales_1.3.0     markdown_1.13    knitr_1.49       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.4.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] labeling_0.4.3    rmarkdown_2.29    compiler_4.4.2

Bibliografia

Johnson, A. A., Ott, M., & Dogucu, M. (2022). Bayes Rules! An Introduction to Bayesian Modeling with R. CRC Press.
Schervish, M. J., & DeGroot, M. H. (2014). Probability and statistics (Vol. 563). Pearson Education London, UK: