5.9 Approssimazione di Laplace

Concludiamo questo capitolo discutendo un problema leggermente diverso da quello della sezione precedente, ma in un certo senso collegato: data una variabile \(X \in \R^d\) di cui sappiamo che la densità non è gaussiana, in quale senso è possibile comunque approssimarla con una gaussiana? In questo modo ad esempio gli strumenti sviluppati per le variabili gaussiane si potrebbero applicare, tenendo conto dell’errore di approssimazione, anche ad altre densità.

Una soluzione particolarmente semplice e spesso efficace è l’approssimazione di Laplace, che consiste nello sviluppare al secondo ordine il logaritmo della densità di \(X\) \[ x \mapsto \log( p(X=x)),\] (supponendo che \(X\) ammetta una densità continua e abbastanza regolare), in un punto di massimo \(x_{{\max}}\) (ossia una moda della densità di \(X\)). Poiché il gradiente si annulla, si trova lo sviluppo \[ \begin{split} & \log( p(X = x) ) \\ & = \log( p(X= x_{{\max}})) + \frac 12 (x-x_{{\max}}) \cdot H(x_{{\max}}) (x-x_{{\max}}) + O( |x-x_{{\max}}|^3), \end{split}\] dove \[ H(x) = \bra{ \frac{ \partial^2 \log (p(X=x))}{\partial x_i \partial x_j}}_{i,j=1}^d\] la matrice delle derivate seconde (detta anche hessiana). Essendo \(x_{{\max}}\) punto di massimo, \(H\) è una matrice (semi-)definita negativa. Supponendo che sia negativa, allora l’approssimazione di Laplace è data dalla densità gaussiana di valor medio \(x_{{\max}}\) e matrice delle covarianze \(- (H(x_{\max}))^{-1}\). L’utilità di questo metodo è che spesso, nel calcolare \(x_{{\max}}\) tramite opportuni metodi numerici si calcola anche la matrice hessiana (ad esempio con il metodo di Newton).

Va evidenziato tuttavia che non è garantita in generale che l’approssimazione sia vicina alla densità di \(X\), in particolare per valori \(x\) lontani da \(x_{{\max}}\). Vediamo degli esempi.

# Iniziamo con una densità non
# gaussiana ma piuttosto simile,

deltax <- 0.01
x <- seq(0, 1, by = deltax)

densita <- x^4 * (1 - x)^4
densita <- densita/sum(densita * deltax)

plot(x, densita, type = "l", lwd = 3, col = miei_colori[1],
  ylab = "densità")

# calcoliamo il minimo dell'opposto del
# logaritmo della densità con la
# funzione nlm() -- questo si può anche
# fare a mano, ma poi lo possiamo
# applicare a casi generali. Dobbiamo
# comunque specificare la funzione
# fuori dall'intervallo [0,1],
# ponendola ad esempio infinita.

log_dens <- function(x) {
  if (x < 0 | x > 1)
    Inf else -log(x^4 * (1 - x)^4)
}

moda <- nlm(log_dens, p = 2/3, hessian = TRUE)

# ricaviamo il punto di massimo e la
# matrice hessiana (qui la derivata
# seconda)

moda$estimate

## [1] 0.4999995

moda$hessian

##      [,1]
## [1,]   32

# aggiungiamo ora al plot la densità
# ottenuta con l'approssimazione di
# Laplace

plot(x, densita, type = "l", lwd = 3, col = miei_colori[1],
  ylab = "densità")
lines(x, dnorm(x, mean = moda$estimate, sd = sqrt(1/moda$hessian)),
  lwd = 3, col = miei_colori[2])


legend("topright", legend = c("originale",
  "approssimata"), fill = miei_colori[1:2])

Basta tuttavia modificare di poco l’esempio sopra per ottenere una cattiva approssimazione, dovuta all’asimmetria della densità.

densita <- x^3 * (1 - x)^5
densita <- densita/sum(densita * deltax)


log_dens <- function(x) {
  if (x < 0 | x > 1)
    Inf else -log(x^3 * (1 - x)^5)
}

moda <- nlm(log_dens, p = 1/2, hessian = TRUE)



plot(x, densita, type = "l", lwd = 3, col = miei_colori[1],
  ylab = "densità")
lines(x, dnorm(x, mean = moda$estimate, sd = sqrt(1/moda$hessian)),
  lwd = 3, col = miei_colori[2])


legend("topright", legend = c("originale",
  "approssimata"), fill = miei_colori[1:2])

Nel caso di densità con più di un punto di massimo locale l’approssimazione è ben peggiore, come è naturale aspettarsi.

x <- seq(-2, 2, by = deltax)
densita <- exp(-(1 - x)^2 * (1 + x)^2)
densita <- densita/sum(densita * deltax)


log_dens <- function(x) {
  if (x < -2 | x > 2) {
    Inf
  } else (1 - x)^2 * (1 + x)^2
}

moda <- nlm(log_dens, p = -1.5, hessian = TRUE)



plot(x, densita, type = "l", lwd = 3, col = miei_colori[1],
  ylim = c(0, 1.2), ylab = "densità")
lines(x, dnorm(x, mean = moda$estimate, sd = sqrt(1/moda$hessian)),
  lwd = 3, col = miei_colori[2])


legend("topright", legend = c("originale",
  "approssimata"), fill = miei_colori[1:2])

Consideriamo infine un esempio nel caso vettoriale. Approssimiamo variabile a valori in \(\R^2\) con densità \[ p( (X,Y ) = (x,y)) \propto (5-(y-\sin(2\pi x))^2 )(1-x^2)(1-y^2)\] per \(x\), \(y \in [-1,1]\) e nulla fuori.

deltax <- 0.05
deltay <- 0.05

x <- seq(-1, 1, by = deltax)
y <- seq(-1, 1, by = deltay)

N_x <- length(x)
N_y <- length(y)


# creiamo una matrice con i valori
# della densità, inizialmente tutti
# nulli

densita <- matrix(nrow = N_x, ncol = N_y)

# definiamo la funzione che calcola la
# densità

densita_funzione <- function(v) {
  (5 - (v[2] - sin(2 * pi * v[1]))^2) *
    (1 - v[1]^2) * (1 - v[2]^2)
}

for (i in 1:N_x) {
  for (j in 1:N_y) {
    densita[i, j] <- densita_funzione(c(x[i],
      y[j]))
  }
}

densita <- densita/sum(densita * deltax *
  deltay)



filled.contour(x, y, densita, color.palette = viridis,
  xlab = "x", ylab = "y")

Figura 5.20: heatmap della densità sopra

Vediamone ora l’approssimazione di Laplace.

log_dens <- function(x) {
  if (abs(x[1]) > 1 | abs(x[2]) > 1) {
    Inf
  } else -log(densita_funzione(x))
}

moda <- nlm(log_dens, p = c(0.5, 0.5), hessian = TRUE)


library(mvtnorm)

m <- moda$estimate

# usiamo la funzione solve per
# calcolare l'inversa della matrice
# hessiana
K <- solve(moda$hessian)


for (i in 1:N_x) {
  for (j in 1:N_y) {
    densita[i, j] <- dmvt(c(x[i], y[j]),
      m, K)
  }
}



filled.contour(x, y, densita, color.palette = viridis,
  xlab = "x", ylab = "y")

Figura 5.21: heatmap dell’approssimazione di Laplace.

5.9.1 Esercizi

Esercizio 5.14 Scrivere l’approssimazione di Laplace per una densità proporzionale a \(x^3 (1-x)^2\) per \(x \in [0,1]\), e nulla altrimenti. Procedere sia analiticamente che tramite opportuni comandi R.