Skip to Tutorial Content

Laplace-Approximation

Die Berechnung des Integrals ist zum Beispiel möglich über Laplace-Integration (die wieder auf Pierre-Simon Laplace zurückgeht). Diese kann allgemein angewandt werden auf Integrale der Form

\[ \int_a^b g(x) = \int_a^b \exp(h(x)) dx \]

wobei \(h(x)\) eine zweimal differenzierbare Funktion ist.

Idee

  • Die Grund-Idee der Laplace-Approximation ist, die Funktion \(g(x)\) als Normalverteilungs-Dichte zu approximieren.
  • Mit \(h(x)=\log(g(x))\) ist \(h(x)\) dann eine quadratische Funktion in \(x\).
  • Wir benutzen die Taylor-Entwicklung bis zum zweiten Term für die Funktion \(h(x)\) um den Punkt \(x_0\):

\[ h(x)\approx h(x_0) + h'(x-x_0)(x-x_0) + \frac{h''(x_0)}{2}(x-x_0)^2 \]

Für die Wahl von \(x_0\) überlegen wir:

  • Wir wollen eigentlich das Integral approximieren.
  • Die Approximation muss also dort gut seinm, wo \(g(x)\) groß ist.
  • Wir wählen also \(x_0=\text{argmax}(g(x))\), das Maximum von \(g(x)\).
  • Das Maximum von \(g(x)\) ist gleichzeitig das Maximum von \(h(x)\).
  • Am Maximum von \(h(x)\) gilt \(h'(x_0)=0\).

Damit erhalten wir:

\[ h(x)\approx h(x_0) + \frac{h''(x_0)}{2}(x-x_0)^2 \]

bzw.:

\[ \begin{aligned} \int g(x) \approx & \int\exp\left(h(x)\right) dx\\ \approx & \int \exp\left(h(x_0) + \frac{h''(x_0)}{2}(x-x_0)^2\right)dx\\ \approx & \,\exp(h(x_0)) \cdot \int \exp\left( \frac{h''(x_0)}{2}(x-x_0)^2\right)dx\\ \approx & \,\exp(h(x_0)) \cdot\sqrt{2\pi(-h''(x_0))}\cdot \\ & \int \frac{1}{\sqrt{2\pi(-h''(x_0))}} \exp\left( \frac{h''(x_0)}{2}(x-x_0)^2\right)dx \end{aligned} \]

Die Funktion im Integral ist nun die Dichte einer normalverteilung mit Erwartungswert \(x_0\) und Varianz \(-h''(x_0)\). Das Integral über die Dichte ist 1.

Damit bleibt als Laplace-Approximation

\[ \int g(x) dx \approx \exp(h(x_0)) \sqrt{2\pi(-h''(x_0))} \]

Algorithmus

Wir brauchen also folgende Schritte:

  1. Finde den Posteriori-Modus
  2. Berechne die zweite Ableitung der Log-Dichte
  3. Setze in die Taylor-Approximation ein

Dabei braucht man für die numerische Optimierung in der Regel sowie die zweite Ableitung.

R

Die folgende R-Funktion berechnet den Logarithmus des Integrals mittels der Laplace-Approximation:

# logpost: Funktion, die Log-Posteriori berechnet
# mode: Startwert
laplace <- function (logpost, mode, y) 
{
  # optim sucht Posteriori-Modus von logpost, gibt zudem zweite Ableitung (hessian) zurück
  fit = optim(mode, logpost, method="Brent", 
    y=y, hessian = TRUE, control = list(fnscale = -1),
    lower = 10^{-6}, upper=10^6)
  mode = fit$par # Modus
  h = -1/fit$hessian # zweite Ableitung am Modus
  int = 1/2 * log(2 * pi) + 0.5 * log(h) + logpost(mode, y) # Taylor-Approximation bis zum zweiten Term
  return(int)
}
  • Es ist meist numerisch günstiger mit der Log-Dichte zu rechnen als mit der richtigen Dichte.

Anwendung

In unserem Beispiel war das erste Modell \(y_t\sim Po(\theta)\) mit Gamma-Priori. Wir brauchen hier also Log-Dichte der Poisson-Verteilung und Log-Dichte der Gamma-Verteilung:

model.1 <- function(theta, y){
  sum(dpois(y, theta, log=TRUE)) + 
    dgamma(theta, shape=280, rate=4, log=TRUE)
}
(log.marg.1 <- laplace(model.1, 70, fire.counts))
##           [,1]
## [1,] -137.7395

Als Startwert für die Suche des Maximums nehmen wir hier den Priori-Erwartungswert 70.

model.2 <- function(theta, y){
  sum(dnorm(y, theta, 6, log=TRUE)) + 
    dgamma(theta, shape=280, rate=4, log=TRUE)
}
model.3 <- function(theta, y){
  sum(dnorm(y, theta, 12, log=TRUE)) + 
    dgamma(theta, shape=280, rate=4, log=TRUE)
}
(log.marg.2 <- laplace(model.2, 70, fire.counts))
##           [,1]
## [1,] -151.0382
(log.marg.3 <- laplace(model.3, 70, fire.counts))
##           [,1]
## [1,] -138.8251

Zurück

Zurück

Laplace-Approximation