8.5 Cenni ai metodi Monte Carlo

Una delle applicazioni principali dei teoremi limite è di utilizzarli insieme alle tecniche di generazione di numeri pseudo-casuali mediante opportuni algoritmi (che non descriviamo nel dettaglio). Tramite semplici comandi in R ad esempio, è possibile simulare variabili aleatorie con densità comuni (uniformi, binomiali, poisson, esponenziali, gaussiane ecc.).

# Il comando per simulare una o più
# variabili aleatorie con una data
# densità si ottiene con il prefisso
# `r` seguito dall'abbreviazione della
# densità. Vediamo ad esempio con le
# gaussiane.

gaussiane <- rnorm(200, mean = 2, sd = 3)

# il comando genera i dati associati
# all'osservazione di 100 variabili
# gaussiane indipendenti.

par(mfrow = c(1, 2))

hist(gaussiane, freq = FALSE, col = miei_colori[1],
  xlab = "valore", ylab = "densità", main = "")

intervallo <- seq(-10, 10, by = 0.01)

lines(intervallo, dnorm(intervallo, mean = 2,
  sd = 3), col = miei_colori[2], lwd = 3)

qqnorm(gaussiane, col = miei_colori[1], pch = 16)
qqline(gaussiane, col = miei_colori[2], lwd = 3)

Figura 8.4: istogramma e qqplot di 200 variabili gaussiane indipendenti

Con opportune estensioni (ad esempio le librerie JAGS o Stan) si può realizzare variabili con densità più complesse. Le tecniche principali consistono nel simulare opportune catene di Markov.

Oltre a fornire dati aleatori ma in un certo senso “controllati” su cui testare esempi e metodi, la simulazione di variabili aleatorie può essere inserita in opportuni algoritmi, che portano il nome di metodi Monte Carlo, in cui anche problemi che non riguardano la probabilità sono affrontabili mediante simulazioni di un gran numero di variabili.

Per fare un esempio, si supponga di dover calcolare l’integrale \[ \int_{\R^d} g(x) p(X=x) d x,\] dove \(g\) e \(p\) sono note e regolari, ma la dimensione \(d\) è molto elevata. Ricorrendo ad una simulazione (purche non abbia costi computazionali onerosi) di \(n\) variabili indipendenti, tutte con densità \(p(X=\cdot)\), allora la legge dei grandi numeri assicura che per \(n\) grande, \[ \int_{\R^d} g(x) p(X=x) d x \approx \frac 1 n \sum_{i=1}^n g(X_i),\] con alta probabilità. Il teorema limite centrale garantisce invece che le oscillazioni \[ \int_{\R^d} g(x) p(X=x) d x- \frac 1 n \sum_{i=1}^n g(X_i)\] saranno circa gaussiane ma con una deviazione standard proporzionale a \(1/\sqrt{n}\). Il vantaggio di questo metodo, rispetto ad esempio ad una integrazione numerica deterministica, è che in dimensione alta (ossia se \(d\) è grande) non è particolarmente più oneroso (mentre in dimensione bassa è possibile fare di meglio).

n <- 1000

gaussiane <- rnorm(n)

g <- gaussiane^2

monte_carlo <- cumsum(g)/1:n

plot(1:n, monte_carlo, type = "l", lwd = 3,
  col = miei_colori[2], ylab = "approssimazione",
  xlab = "numero simulazioni")

abline(h = 1, col = miei_colori[1], lwd = 3)

Figura 8.5: Calcolo del momento secondo di una gaussiana standard tramite media empirica di un campione indipendente.