Julia Community 🟣: Martin Roa Villescas

Naïve k-means

Martin Roa Villescas — Tue, 28 Jun 2022 19:16:26 +0000

This is the first of a series of posts on k-means clustering. In this post, we will focus on k-means clustering in its simplest form, commonly referred to as naïve k-means. In future posts, we will gradually improve this version by better understanding its inner workings and looking at it from a different angle.

The goal of k-means clustering is to partition a set of $N$ observations ${x1,…,xN}\{ \mathbf{x}_1, \ldots, \mathbf{x}_N \}$ into $K$ clusters. An observation $xn\mathbf{x}_n$ corresponds to a point in a multidimensional space. Intuitively, each cluster $k$ can be thought of as a group of points that lie close to each other around its mean $μk.\bm{\mu}_k.$ The goal of k-means clustering is to find an assignment of points to clusters ${ k_n \}$ that minimizes the sum of the squared Euclidean distances of each point $xn\mathbf{x}_n$ to its closest mean $μk.\bm{\mu}_k.$ We can do this through an iterative procedure involving two steps:

Assign each point $xn\mathbf{x}_n$ to its nearest cluster $k_n$ with mean $μk\bm{\mu}_k$ : $mink(μk−xn)2k_n = \text{arg}~\text{min}_k (\bm{\mu}_k - \mathbf{x_n})^2$
Update the cluster means ${μk}\{ \bm{\mu}_k \}$ :

\bm{\mu_k} = \frac{\sum\limits_{\bm{x_n} \in k} \bm{x_n}}{\sum\limits_{\bm{x_n} \in k} 1}

This two-step iterative procedure is repeated until convergence. Note that before we can start iterating, we need to choose initial values for the means $μk.\bm{\mu}_k.$ One common approach is to set them to a random subset of $K$ observations. For the moment, we suppose that the value of $K$ is given.

Let us now look at an implementation using the Julia language.

We start start by importing the necessary packages: Plots.jl to visualize the result, RDatasets.jl to load the Iris flower data set, Distances.jl to use its squared Euclidean distance implementation and, Statistics.jl to use its mean implementation. We then load two features of the dataset into a matrix. For illustration purposes we will use two features. Note, however, that k-means clustering works with data of any number of dimensions. Finally, we choose an arbitrary value for $K .$

using Plots, RDatasets, Distances, Statistics
pyplot(leg=false, ms=6, border=true, fontfamily="Calibri",
       title="Naïve k-means")             # plot defaults
iris = dataset("datasets", "iris");       # load dataset
𝐱ₙ = collect(Matrix(iris[:, 1:2])');      # select features
K = 3                                     # number of clusters

Let's define a function to visualize the cluster points ${xn}\{ \mathbf{x}_n \}$ and the means ${μk}.\{ \bm{\mu}_k \}.$ Let's also represent the assignment of points to clusters ${ k_n \}$ using the marker's color.

function plot_clusters(𝐱ₙ, 𝛍ₖ, kₙ)
  plt = scatter(𝐱ₙ[1,:], 𝐱ₙ[2,:], c=kₙ)
  scatter!(𝛍ₖ[1,:], 𝛍ₖ[2,:], m=(:xcross,10), c=1:size(𝛍ₖ,2))
  return plt
end

Last but not least, an elegant but naïve k-means implementation:

# Initialize the means for each cluster
𝛍ₖ = 𝐱ₙ[:, rand(axes(𝐱ₙ, 2), K)]

anim = @animate for _ in 1:15

  # Find the squared distance from each sample to each mean
  d = pairwise(SqEuclidean(), 𝐱ₙ, 𝛍ₖ, dims=2)

  # 1. Assign each point to the cluster with the nearest mean
  kₙ = findmin(d, dims=2)[2] |> x -> map(i -> i.I[2], x)

  # 2. Update the cluster means
  for k in 1:K
    𝐱ₙₖ = 𝐱ₙ[:, dropdims((kₙ .== k), dims=2)]
    𝛍ₖ[:, k] = mean(𝐱ₙₖ, dims=2)
  end

  plot_clusters(𝐱ₙ, 𝛍ₖ, kₙ)

end

gif(anim, joinpath(@__DIR__, "anim.gif"), fps=2)

K-means clustering is one of the simplest and most popular unsupervised machine learning algorithms used to discover underlying patterns by grouping similar data. In this post, we studied the mechanics of this method in its simplest form. However, a number of questions remain unanswered. For example, how did we arrive to the equations describing the two-step iterative procedure? Is it guaranteed to converge? What are the shortcomings of k-means clustering as implemented here? Can we give a meaningful interpretation to this problem that gives us a deeper understanding as to what is actually happening in the iterative procedure? If so, can such interpretation provide us with tools to modify and improve the mechanism presented here? Stay tuned.

The product of two Gaussian PDFs is Gaussian. Is it?

Martin Roa Villescas — Wed, 22 Jun 2022 17:36:03 +0000

A few years back, I followed a course on Bayesian machine learning. Very early in the course, we were introduced to the Gaussian probability density function (PDF). We were presented with the equation that describes the univariate case. Soon after, we arrived to the following conclusion:

The multiplication of two Gaussian PDFs yields another (unnormalized) Gaussian.

This claim is quite relevant for Bayesian machine learning given that the product of two distributions is exactly what's going on in the numerator of Bayes' theorem:

\overbrace{p(\lambda \mid D)}^{\text{posterior}} = \frac{\overbrace{p(D \mid \lambda)}^{\text{likelihood}} \, \overbrace{p(\lambda)}^{\text{prior}}}{\underbrace{p(D)}_{\text{evidence}}},

The professor even walked us through the proof (e.g. like in here).

However, in my head, there was something that didn't click. I remember I asked myself 'Wait, you mean that if we multiply the following two Gaussian PDFs, which are quite displaced with respect each other, point by point, we still get another bell-shaped curve as a result? Don't they cancel out each other? Yes, they must cancel out...'

So yeah. I was skeptical. It was basically a battle between my intuition and mathematics.

So in order to see who was right, I figured that plotting the result of this operation would tell. In the next figure, the top and middle rows correspond to the two operand Gaussians. The bottom row contains the 2 operand Gaussians and their product.

'I knew it,' I said when I first saw this result. If you look at the bottom plot closely, you can see a green line along the x axis, which is the result of the product. 'It cancels out!'.

Soon after I thought 'Hmmm, wait. What if the result is actually Gaussian but too small to see. I need a way to zoom in.' Luckily, Plots.jl makes it very easy to animate this kind of things. So I did.

using Plots, LaTeXStrings, Printf, Distributions
pyplot(f=true, fa=0.3, grid=false, xticks=false, size=(600, 400),
       legend=false, xlabel=L"\theta", ylabel=L"P(\theta)",
       fontfamily="CMU Serif", yformatter=y->@sprintf("%.1E",y))

θ = range(-10, 10, length=500)

d1 = pdf.(Normal(-5, 1), θ)
d2 = pdf.(Normal(5, sqrt(2)), θ)

p1 = plot(d1, yticks=[maximum(d1)])      # top
p2 = plot(d2, yticks=[maximum(d2)], c=2) # middle

ymaxr = [0.4/i^4.5 for i=1:50] |> x->vcat(x,repeat([last(x)],20))

anim = @animate for ymax in ymaxr
  p3 = plot(d1)                          # bottom
  plot!(d2)
  plot!(d1 .* d2, yticks=[ymax*0.97], ylims=(0,ymax))
  plot(p1, p2, p3, layout=(3,1))         # all
end

gif(anim, joinpath(@__DIR__, "anim.gif"), fps=10)

So yeah,

Mathematics: 1
Martin: 0

The Laplace approximation (part 2)

Martin Roa Villescas — Fri, 17 Jun 2022 19:43:04 +0000

This is the second post in a two-part series on the Laplace transformation. In part 1, we saw how a well-behaved general function $f (x)$ can be approximated with an unnormalized normal distribution. In this post, I will extend this idea by showing how the Laplace approximation is used in the context of Bayesian inference. I present an example implemented in Julia that compares the Laplace approximation of a probability distribution with its exact counterpart.

Bayesian inference

Bayesian inference is a method used to update the probability for a hypothesis as more information becomes available. At the core of this process lies Bayes' theorem, a mathematical formula that transforms the prior distribution of a variable into a posterior distribution in light of new or additional data:

\overbrace{p(\lambda \mid D)}^{\text{posterior}} = \frac{\overbrace{p(D \mid \lambda)}^{\text{likelihood}} \, \overbrace{p(\lambda)}^{\text{prior}}}{\underbrace{p(D)}_{\text{evidence}}},

where $λ\lambda$ is the variable of interest and $D$ the set of observations. Using the chain rule and the law of total probability we can rewrite Bayes' theorem as

\overbrace{p(\lambda \mid D)}^{\text{posterior}} = \frac{\overbrace{p(D, \lambda)}^{\text{joint likelihood}} }{\underbrace{\int_{\lambda}p(D, \lambda) \, \mathrm{d}\lambda}_{\text{evidence}}}.

In many practical applications, the integral in the equation above is intractable, i.e., impossible to solve analytically. Therefore, we have to resort to approximation methods to compute it. The Laplace approximation is one example.

The Laplace approximation

In contrast with part 1, where the goal was to find the Laplace approximation $q (x)$ of a general function $f (x)$ , here we want to approximate a probability distribution, which by definition integrates to one over its entire domain. This means that in addition to finding the unnormalized distribution $q (x)$ , we also need to calculate its normalizing constant $C$ . I.e. our goal is to find

p_{\mathcal{L}}(\lambda \mid D) = \frac{1}{C} f(\lambda),

which is nothing more than a disguised form of Bayes' theorem:

p_{\mathcal{L}}(\lambda \mid D) = \frac{p_{\mathcal{L}}(D, \lambda)}{p_{\mathcal{L}}(D)}.

In the following sections, we will derivate the equations for $pL(D,λ)p_{\mathcal{L}}(D, \lambda)$ and $pL(D)p_{\mathcal{L}}(D)$ according to the Laplace approximation.

Even though the derivation of these equations might seem tedious, the conclusion is succinct and elegant, and is all you need in practice.

Step 1: approximating the joint likelihood $pL(D,λ)p_{\mathcal{L}}(D, \lambda)$

We start by finding the mode $λ0\lambda_0$ of the joint likelihood $\lambda)$ , or equivalently, of the joint log-likelihood $ln⁡p(D,λ)\ln p(D, \lambda)$ , which can be done either analytically or numerically. The latter option is demonstrated later in the example.

Then, we make a second-order Taylor expansion of $\lambda_0)$ around the mode $λ0\lambda_0$ , like we did in part 1. The procedure is repeated here for completeness.

To simplify the notation, let $f(λ)=pL(D,λ)f(\lambda) = p_{\mathcal{L}}(D, \lambda)$ and $h(λ)=ln⁡f(λ)h(\lambda) = \ln f(\lambda)$ . Consider the second-order Taylor expansion of $h(λ)h(\lambda)$ around $λ0\lambda_0$ :

h(\lambda) \approx h(\lambda_0) + h'(\lambda_0)(\lambda-\lambda_0) + \frac{1}{2}h''(\lambda_0)(\lambda-\lambda_0)^2

Note that $h′(λ0)=dln⁡f(λ)dλ∣λ=λ0=f′(λ0)f(λ0)=0h'(\lambda_0) = \frac{\mathrm{d}\ln f(\lambda)}{\mathrm{d}\lambda}\Bigr|_{\lambda = \lambda_0} = \frac{f'(\lambda_0)}{f(\lambda_0)} = 0$ , since $f′(λ0)=0f'(\lambda_0) =0$ , so we have:

h(\lambda) \approx h(\lambda_0) + \frac{1}{2}h''(\lambda_0)(\lambda-\lambda_0)^2.

Applying the exponential function on both sides yields:

\begin{aligned} \exp \left( h(\lambda) \right) &\approx \exp\left( h(\lambda_0) + \frac{1}{2}h''(\lambda_0)(\lambda-\lambda_0)^2 \right) \\ \exp \left( \ln f(\lambda) \right) &\approx \exp\left( h(\lambda_0) \right) \exp \left( \frac{1}{2}h''(\lambda_0)(\lambda-\lambda_0)^2 \right) \\ f(\lambda) &\approx f(\lambda_0) \exp \left( \frac{1}{2}h''(\lambda_0)(\lambda-\lambda_0)^2 \right), \end{aligned}

Hence,

\begin{aligned} p_{\mathcal{L}}(D, \lambda) &\approx p(D, \lambda_0) \exp \left( \frac{1}{2}h''(\lambda_0)(\lambda-\lambda_0)^2 \right), \end{aligned}

Step 2: approximating the evidence $pL(D)p_{\mathcal{L}}(D)$

Again, to simplify the notation, let $f(λ)=pL(D,λ)f(\lambda) = p_{\mathcal{L}}(D, \lambda)$ and $p_{\mathcal{L}}(D)$ . According to the law of total probability, we have that

\int_{\lambda} f(\lambda) \, \mathrm{d}\lambda.

Substituting the result that we derived above for $f(λ)f(\lambda)$ , we obtain

\begin{aligned} C &\approx f(\lambda_0) \int_{\lambda} \exp \left( \frac{1}{2}h''(\lambda_0)(\lambda-\lambda_0)^2 \right) \, \mathrm{d}\lambda \\ &\approx f(\lambda_0) \begingroup\color{RoyalBlue} \int_{\lambda} \exp \left( -\frac{(-h''(\lambda_0))}{2}(\lambda-\lambda_0)^2 \right) \, \mathrm{d}\lambda \endgroup. \end{aligned}

Recall that the probability density function (PDF) of a normally distributed random variable $X$ with mean $μ\mu$ and variance $σ2\sigma^2$ , denoted as $\sim \mathcal{N}(\mu, \sigma^2)$ , is given by

\frac{1}{\sqrt{2\pi}\sigma}\exp\left(-\frac{1}{2}\frac{(x-\mu)^2}{\sigma^2}\right),

or equivalently by

\frac{\sqrt{\gamma}}{\sqrt{2\pi}} \exp \left(-\frac{\gamma}{2}(x-\mu)^2\right),

where the precision $γ=1σ2\gamma= \frac{1}{\sigma^2}$ .

Integrating both sides with respect to $x$ yields

\begin{aligned} \int_x p(X) \, \mathrm{d}x &= \frac{\sqrt{\gamma}}{\sqrt{2\pi}} \int_x \exp \left(-\frac{\gamma}{2}(x-\mu)^2\right) \, \mathrm{d}x \\ 1 &= \frac{\sqrt{\gamma}}{\sqrt{2\pi}} \int_x \exp \left(-\frac{\gamma}{2}(x-\mu)^2\right) \, \mathrm{d}x \\ \frac{\sqrt{2\pi}}{\sqrt{\gamma}} &= \begingroup\color{RoyalBlue} \int_x \exp \left(-\frac{\gamma}{2}(x-\mu)^2\right) \, \mathrm{d}x \endgroup. \end{aligned}

If we compare the two equations in blue above, we can conclude that

\int_{\lambda} \exp \left( -\frac{(-h''(\lambda_0))}{2}(\lambda-\lambda_0)^2 \right) \, \mathrm{d}\lambda = \frac{\sqrt{2\pi}}{\sqrt{-h''(\lambda_0)}}.

Hence,

\begin{aligned} C &\approx f(\lambda_0) \frac{\sqrt{2\pi}}{\sqrt{-h''(\lambda_0)}}, \\ p_{\mathcal{L}}(D) &\approx p(D, \lambda_0) \frac{\sqrt{2\pi}}{\sqrt{-h''(\lambda_0)}}. \end{aligned}

Step 3: bringing it all together

Recall that we seek to approximate the posterior distribution $p(λ∣D)p(\lambda \mid D)$ with its Laplace approximation $pL(λ∣D)=pL(D,λ)pL(D)p_{\mathcal{L}}(\lambda \mid D) = \frac{p_{\mathcal{L}}(D, \lambda)}{p_{\mathcal{L}}(D)}$ .

Using a second-order Taylor expansion, we derived that

p_{\mathcal{L}}(D, \lambda) \approx p(D, \lambda_0) \exp \left( \frac{1}{2}h''(\lambda_0)(\lambda-\lambda_0)^2 \right),

We also showed that the normalizing constant is

p_{\mathcal{L}}(D) \approx p(D, \lambda_0) \frac{\sqrt{2\pi}}{\sqrt{-h''(\lambda_0)}}.

Substituting these two results into $pL(λ∣D)=pL(D,λ)pL(D)p_{\mathcal{L}}(\lambda \mid D) = \frac{p_{\mathcal{L}}(D, \lambda)}{p_{\mathcal{L}}(D)}$ and rearranging the equation, we obtain

\begin{aligned} p_\mathcal{L}(\lambda \mid D) &= \frac{\sqrt{-h''(\lambda_0)}}{p(D, \lambda_0)\sqrt{2\pi}} p(D, \lambda_0) \exp \left( \frac{1}{2} h''(\lambda_0)(\lambda - \lambda_0)^2 \right) \\ &= \frac{1}{\sqrt{2\pi}\sqrt{(-h''(\lambda_0))^{-1}}} \exp \left( -\frac{1}{2} \frac{(\lambda - \lambda_0)^2}{(-h''(\lambda_0))^{-1}} \right), \end{aligned}

or more succinctly

\boxed{p_\mathcal{L}(\lambda \mid D) = \mathcal{N} \left( \lambda_0, (-h''(\lambda_0))^{-1} \right)} \, ,

where

\boxed{h''(\lambda_0) = \frac{\mathrm{d}^2 \ln p(D, \lambda)}{\mathrm{d}\lambda^2}\Bigr|_{\lambda = \lambda_0}} \, .

In conclusion, in the context of Bayesian inference, the Laplace approximation of a posterior distribution is a normal distribution centered around the mode of the joint likelihood and has a precision equal to minus the second derivative of the joint log-likelihood evaluated at its mode.

Example: Poisson data with a gamma prior

Let's look at a simple example of a model that has an exact solution. This will allow us to compare the Laplace approximation of a posterior distribution with its exact counterpart. Our model assumes that the data is drawn from a Poisson distribution with a gamma distributed rate parameter $λ\lambda$ . I.e. the model is defined as

\begin{aligned} Y &\sim \text{Poisson}(\lambda) \\ \lambda &\sim \text{Gamma}(\alpha, \theta). \end{aligned}

The objective is to calculate $p(λ∣Y)p(\lambda \mid Y)$ . The exact solution is given by

p(\lambda \mid \{y_i\}) \sim \text{Gamma}\left(\alpha + \sum_{i=1}^n y_i, \frac{\theta}{n\theta + 1} \right),

where ${y_i\}$ are the set of observed values.

We start by importing the necessary packages and by setting some default options for the plotting package:

using Plots, LaTeXStrings         # used to visualize the result
using SpecialFunctions: gamma     # used to define the gamma PDF
using Optim: maximize, maximizer  # used to find x₀
using ForwardDiff: derivative     # used to find h''(x₀)
default(fill=true, fillalpha=0.3, xlabel="λ") # plot defaults

Next, we define the normal, Poisson, and gamma PDFs:

Normal(μ, σ²) = x -> (1/(sqrt(2π*σ²))) * exp(-(x-μ)^2/(2*σ²))
Poisson(Y::Int) = λ -> (λ^Y) / factorial(Y) * exp(-λ)
Gamma(α, θ) = λ -> (λ^(α-1) * exp(-λ/θ)) / (θ^(α) * gamma(α))

The next step is to specify the model. To do so, we need to define the prior $p(λ)p(\lambda)$ and the likelihood $\mid \lambda)$ according to the specification above. To define the prior, we need to set the hyperparameters of the gamma distribution $α\alpha$ and $θ\theta$ . We will set them both to 3. To specify the likelihood, we use a Poisson distribution with a placeholder $Y$ that we can use to enter observed values. The model is now fully specified. However, let's define the joint log-likelihood $ln⁡p(D,λ)\ln p(D, \lambda)$ , (the log of the product of the likelihood and the prior) which is a quantity needed by the Laplace approximation:

prior = Gamma(3, 3)
likelihood(Y) = Poisson(Y)
joint_log_likelihood(Y) = λ -> log(likelihood(Y)(λ) * prior(λ))

Suppose we observe $y = 2$ . The exact posterior is given by:

y = 2
exact_posterior = Gamma(prior.α + y, prior.θ / (prior.θ + 1))

We can visualize the prior and exact posterior:

λ = range(0, 12, length=500)
plot(λ, prior, lab=L"p(\lambda)")
plot!(λ, exact_posterior, lab=L"p(\lambda|D)")

Let's now find the Laplace approximation of the exact posterior using the results we derived. The first step is to find the mode $λ0\lambda_0$ of the joint log-likelihood $ln⁡p(Y,λ)\ln p(Y, \lambda)$ . We can do so using the maximize function from the package Optim.jl. The arguments expected by maximize are the function to be maximized and the endpoints of the interval along which the optimization will be performed. The returned object is a structure with several information from which we can extract the mode using the maximizer function.

λ₀ = maximize(joint_log_likelihood(y), first(λ), last(λ)) |> maximizer

The final step is to find $ln⁡h′′(λ0)\ln h''(\lambda_0)$ , i.e., the second derivative of the joint log-likelihood $ln⁡p(Y,λ)\ln p(Y, \lambda)$ at its mode $λ0\lambda_0$ . We use the derivative function from the ForwardDiff.jl package:

hꜛꜛ(λ) = derivative(λ -> derivative(joint_log_likelihood(y), λ), λ)

Finally, we bring together the pieces to define the Laplace approximation:

laplace_approx = Normal(λ₀, -(hꜛꜛ(λ₀))^(-1))

We can now compare the Laplace approximation to the exact posterior:

plot!(λ, laplace_approx, lab=L"p_{\mathcal{L}}(\lambda|D)", title=L"\textrm{Laplace~approximation}")

A few remarks to conclude this post. The Laplace approximation replaces the problem of integration with that of optimization. This is important because intractable integrals are common in practical applications. On the other hand, even when the integral can be solved, optimization is often faster to compute. However, the Laplace approximation has several limitations: if the posterior distribution being approximated is not unimodal, or does not have it's mass concentrated around the mode, or is not twice-differentiable, then the result will be inaccurate. This is worrying because in many applications we don't have access to the exact posterior, which means that we don't have a reference to test the accuracy of the approximation. In any case, the Laplace approximation has been proven to be effective not only in the context of Bayesian inference, but also in many other mathematical domains that exhibit hard-to-compute integrals.

Credit to @OswaldoGressani for his awesome YouTube video

The Laplace approximation (part 1)

Martin Roa Villescas — Sun, 12 Jun 2022 11:09:03 +0000

This is the first of a series of posts that I'm planning to upload during the upcoming weeks. The covered material is based on the notes that I've gathered during the time I've been using Julia for scientific computing.

The general idea behind the Laplace approximation is to approximate a well-behaved uni-modal function with an unnormalized Gaussian distribution. "Well-behaved" means that this function should be twice-differentiable and that most of its area under the curve should be concentrated around the mode.

To approximate $f (x)$ , we first make a second-order Taylor expansion around its mode $x_0$ .

Let $\ln f(x)$ and consider the second-order Taylor expansion of $h (x)$ around $x_0$ :

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

Note that $h′(x0)=dln⁡f(x)dx∣x=x0=f′(x0)f(x0)=0h'(x_0) = \frac{\mathrm{d}\ln f(x)}{\mathrm{d}x}\Bigr|_{x = x_0} = \frac{f'(x_0)}{f(x_0)} = 0$ , since $f'(x_0) =0$ , so we have:

\approx h(x_0) + \frac{1}{2}h''(x_0)(x-x_0)^2.

Applying the exponential function on both sides yields:

\begin{aligned} f(x) &\approx \exp\left( h(x_0) + \frac{1}{2}h''(x_0)(x-x_0)^2 \right) \\ &\approx \exp\left( h(x_0) \right) \exp \left( \frac{1}{2}h''(x_0)(x-x_0)^2 \right) \\ &\approx f(x_0) \exp \left( \frac{1}{2}h''(x_0)(x-x_0)^2 \right), \end{aligned}

which corresponds to a Gaussian distribution up to a scaling constant.

As a result, to approximate $f (x)$ we need two things:

The mode $x0\bm{x_0}$ of either $f (x)$ or $h (x)$ (which is the same).
The second derivative of $h (x)$ at $x_0$ $h′′(x0)\bm{h''(x_0)}$ .

Let's look at an example implemented in Julia.

The task is to find the Laplace approximation $q (x)$ of the function:

f(x) = (4-x^2)\exp(-x^2)

We start by importing the packages that we will be needing:

using Plots, LaTeXStrings         # used to visualize the result
using Optim: maximize, maximizer  # used to find x₀
using ForwardDiff: derivative     # used to find h''(x₀)

We continue by defining the function $f (x)$ :

f(x) = (4-x^2)*exp(-x^2)

The next step is to find the mode $x_0$ of $f (x)$ , i.e., the argument that maximizes the function $f (x)$ . For this, we will use the maximize function from Optim.jl, a Julia optimization package. The arguments expected by maximize are the function to be maximized and the endpoints of the interval along which the optimization will be performed. The returned object is a structure with several information from which we can extract the mode using the maximizer function.

x = range(-2, 2, length=100)
x₀ = maximize(x -> f(x), first(x), last(x)) |> maximizer

The final step is to find $h''(x_0)$ , i.e., the second derivative of $h (x)$ at its mode $x_0$ . ForwardDiff.jl is performant Julia package that implements the derivative method to perform automatic differentiation:

h(x) = log(f(x))
hꜛꜛ(x) = derivative(x -> derivative(h, x), x) # second derivative

Now we can bring everything together to define the Laplace approximation $q (x)$ of our original function $f (x)$ :

q(x) = f(x₀) * exp((1/2) * hꜛꜛ(x₀) * (x-x₀)^2)

Note the resemblance of this function with the one we derived above. Julia is pretty neat, huh?

Finally, let's visualize the result:

plot( x, f, lab=L"f(x)= (4-x^2) \cdot e^{-x^2}", xlab=L"x", title=L"\textrm{Laplace~approximation}")
plot!(x, q, lab=L"q(x) = %$(f(x₀)) \cdot e^{\frac{%$(hꜛꜛ(x₀))}{2} x^2}")
sticks!([x₀], [f(x₀)], m=:circle, lab="")

Pretty good approximation for this particular function!

To conclude, we showed how the Laplace approximation is a useful technique to approximate a well-behaved uni-modal function with an unnormalized Gaussian distribution, which is a well-understood and thoroughly studied quantity. Note, however, that if the function being approximated is not unimodal, or does not have it's mass concentrated around the mode, or is not twice-differentiable, then the Laplace approximation is not an adequate option. In part 2, I will go over the role that the Laplace approximation plays in the context of Bayesian inference.