Learning Objectives

By the end of this lesson, you will be able to:

• Understand why approximations are needed in Bayesian statistics.
• Derive the Laplace approximation.
• Approximate difficult posterior distributions.
• Use Laplace approximation in healthcare prediction problems.
• Implement Laplace approximation in Python.

Why Do We Need Approximations?

Suppose we observe patient data:

and want the posterior distribution

Using Bayes’ theorem:

The denominator

is often impossible to compute analytically.

Laplace approximation replaces the posterior with a Gaussian distribution near its peak.

Core Idea

Suppose

where

Let

be the posterior mode.

Use a second-order Taylor expansion:

because

at the maximum.

Because $(\hat{\theta})$ is the value where $(h(\theta))$ is maximum.

The Taylor expansion around $(\hat{\theta})$ is:

But at the maximum:

So the middle term disappears:

Since it is a maximum:

Now exponentiate:

So:

This has the same shape as a Normal distribution:

Match the terms:

Therefore:

Then

This is the Laplace approximation.

Visual Interpretation

The true posterior might be irregular.

Near its peak, we fit a parabola.

Since the exponential of a parabola is a Gaussian,

we obtain a normal approximation.

Healthcare Example

Suppose:

• $y_i=1$ if patient survives surgery.
• $y_i=0$ otherwise.

Logistic model:

Prior:

The posterior has no closed form.

Laplace approximation provides:

allowing confidence intervals and prediction.

Pharmaceutical Example

A drug company studies treatment response.

Parameter:

Posterior:

may be difficult to compute.

Laplace approximation quickly provides:

• Posterior mean
• Posterior variance
• Credible intervals

without expensive simulations.

Connection to GLMs

Laplace approximation is heavily used in:

• Logistic regression
• Poisson regression
• Mixed models
• Bayesian hierarchical models

Many software packages use Laplace internally.

Examples:

• INLA
• glmer
• TMB

Python Example

import numpy as np
from scipy.optimize import minimize
from scipy.stats import norm

y = np.array([1,1,1,0,1,0,1,1,0,1])

def log_posterior(beta):

    p = 1/(1+np.exp(-beta))

    loglik = np.sum(
        y*np.log(p) +
        (1-y)*np.log(1-p)
    )

    logprior = norm.logpdf(
        beta,
        loc=0,
        scale=np.sqrt(10)
    )

    return loglik + logprior

result = minimize(
    lambda b: -log_posterior(b),
    x0=0
)

beta_hat = result.x[0]

eps = 1e-5

hessian = (
    log_posterior(beta_hat+eps)
    -
    2*log_posterior(beta_hat)
    +
    log_posterior(beta_hat-eps)
)/(eps**2)

variance = -1/hessian

print(beta_hat)
print(variance)

import numpy as np
from scipy.optimize import minimize
from scipy.stats import norm

y = np.array([1,1,1,0,1,0,1,1,0,1])

def log_posterior(beta):

    p = 1/(1+np.exp(-beta))

    loglik = np.sum(
        y*np.log(p) +
        (1-y)*np.log(1-p)
    )

    logprior = norm.logpdf(
        beta,
        loc=0,
        scale=np.sqrt(10)
    )

    return loglik + logprior

result = minimize(
    lambda b: -log_posterior(b),
    x0=0
)

beta_hat = result.x[0]

eps = 1e-5

hessian = (
    log_posterior(beta_hat+eps)
    -
    2*log_posterior(beta_hat)
    +
    log_posterior(beta_hat-eps)
)/(eps**2)

variance = -1/hessian

print(beta_hat)
print(variance)

Explaining the code

Step 1: Laplace Approximation Idea

Suppose:$$p(\beta \mid D)\propto p(D\mid \beta )p(\beta )$$

Take logs:$$h(\beta )=\log p(D\mid \beta )+\log p(\beta )$$

This is exactly your:

def log_posterior(beta):

def log_posterior(beta):

function.

Step 2: Find the Peak

Laplace says:

Find the value of $\beta$ that maximizes the log posterior.

That’s:$$\widehat{\beta }=\arg \max h(\beta )$$

which you obtain using:

result = minimize(    lambda b: -log_posterior(b),    x0=0)

result = minimize(    lambda b: -log_posterior(b),    x0=0)

because minimizing the negative is equivalent to maximizing.

Step 3: Taylor Expansion

Near the maximum:$$h(\beta )\approx h(\widehat{\beta })+h^{\mathrm{\prime }}(\widehat{\beta })(\beta -\widehat{\beta })+\frac{1}{2}{h}^{\mathrm{\prime }\mathrm{\prime }}(\widehat{\beta })(\beta -\widehat{\beta }{)}^2$$

At the maximum:$$h^{\mathrm{\prime }}(\widehat{\beta })=0$$

so:$$h(\beta )\approx h(\widehat{\beta })+\frac{1}{2}{h}^{\mathrm{\prime }\mathrm{\prime }}(\widehat{\beta })(\beta -\widehat{\beta }{)}^2$$

This is the key Laplace step.

Step 4: Why We Need the Hessian

Notice that the approximation depends on:$$h^{\mathrm{\prime }\mathrm{\prime }}(\widehat{\beta })$$

the curvature at the peak.

If:$$h^{\mathrm{\prime }\mathrm{\prime }}(\widehat{\beta })=-100$$

the peak is narrow.

If:$$h^{\mathrm{\prime }\mathrm{\prime }}(\widehat{\beta })=-1$$

the peak is wide.

Therefore the curvature determines the variance.

Step 5: Your Code Computes That Curvature

The second derivative formula is:$$h^{\mathrm{\prime }\mathrm{\prime }}(x)\approx \frac{h(x+\epsilon )-2h(x)+h(x-\epsilon )}{{\epsilon }^2}$$

This is exactly:

hessian = (    log_posterior(beta_hat+eps)    -    2*log_posterior(beta_hat)    +    log_posterior(beta_hat-eps))/(eps**2)

hessian = (    log_posterior(beta_hat+eps)    -    2*log_posterior(beta_hat)    +    log_posterior(beta_hat-eps))/(eps**2)

So the code is estimating:$$h^{\mathrm{\prime }\mathrm{\prime }}(\widehat{\beta })$$

numerically.

Step 6: Convert Curvature to Variance

Laplace says:$$p(\beta \mathrm{\mid }D)\approx N(\widehat{\beta },-\frac{1}{h^{\mathrm{\prime }\mathrm{\prime }}(\widehat{\beta })})$$

The variance is:$$V=-\frac{1}{h^{\mathrm{\prime }\mathrm{\prime }}(\widehat{\beta })}$$

which is exactly:

variance = -1/hessian

variance = -1/hessian

Geometric Interpretation

The second derivative measures how sharply it bends at the top.

• Sharp peak → large negative curvature → small variance.
• Flat peak → small negative curvature → large variance.

Laplace simply says:

“Replace the posterior by the normal distribution that has the same peak and curvature.”

The Hessian is the quantity that captures that curvature.

In Multiple Dimensions

Later, when you learn Bayesian logistic regression with many parameters:$$\beta =({\beta }_1,\dots ,{\beta }_p)$$

the second derivative becomes a matrix:$$H=\frac{{\mathrm{\partial }}^{2}h}{\mathrm{\partial }\beta \mathrm{\partial }{\beta }^T}$$

called the Hessian matrix.

Then Laplace becomes:$$p(\beta \mathrm{\mid }D)\approx N(\widehat{\beta },-H^{-1})$$

Your one-dimensional code is the simplest possible version of this idea. The numerical second derivative is estimating the Hessian, and the inverse Hessian gives the posterior variance. That’s the central mathematical connection between the code and the Laplace approximation.

Advantages

• Extremely fast.
• Easy implementation.
• Useful for large datasets.

Limitations

• Assumes posterior is approximately Gaussian.
• Poor for multimodal distributions.
• Poor for highly skewed distributions.

References

• Bishop, Pattern Recognition and Machine Learning.
• Murphy, Machine Learning.
• Gelman et al., Bayesian Data Analysis.
• Rue et al., INLA Methodology.