TruncatedMVN

Julia reimplementation of a truncated multivariate normal distribution with perfect sampling.

Distribution and sampling method by [1].

Code based on MATLAB implementation by Zdravko Botev and Python implementation by Paul Brunzema. These may be found here: MATLAB by Botev, Python by Brunzema

Exports 1 struct + its constructor and 1 function:

• TruncatedMVN.TruncatedMVNormal: Struct and constructor for initializing the distribution.

• TruncatedMVN.sample: Exact sampling from the distribution.

Documentation for these can be found in the Public API section.

Example : Simulation of Gaussian variables subject to linear restrictions

This example is adapted from the MATLAB package image and the R package TruncatedNormal vignette.

Suppose we wish to simulate a bivariate vector $\mathit{X}\sim \mathcal{N}(\mathit{\mu },\mathbf{\Sigma })$ conditional on $X_1-X_2<-6$. Setting

$$\mathbf{A}=\left(\begin{array}{cc}1& -1\\ 0& 1\end{array}\right),$$

we can simulate $\mathit{Y}\sim \mathcal{N}(\mathbf{A}\mathit{\mu },\mathbf{A}\mathbf{\Sigma }{\mathbf{A}}^{\mathrm{\top }})$ conditional on $\mathit{l}\le \mathit{Y}\le \mathit{u}$ and then get $\mathit{X}={\mathbf{A}}^{-1}\mathit{Y}$.

using Random, Distributions
using CairoMakie, LaTeXStrings
using TruncatedMVN

# Problem setup
Random.seed!(1234)
Σ = [1.0 0.9;
     0.9 1.0]
μ = [-3.0, 0.0]
ub = [-6.0, Inf]
lb = [-Inf, -Inf]
A = [1.0 -1.0;
     0.0  1.0]

n = 1000

# Sample from truncated MVN under A*x ≤ ub
tmvn = TruncatedMVNormal(A * μ, A * Σ * A', lb, ub)
Y_truncated = TruncatedMVN.sample(tmvn, n)          # 2×n matrix
X_conditional = A \ Y_truncated                       # 2×n matrix

# Unconstrained samples for comparison
X_unconditional = rand(MvNormal(μ, Σ), n)     # 2×n matrix

# Plot
with_theme(theme_latexfonts(), fontsize = 18) do
    fig = Figure(size=(960, 540))
    ax = Axis(fig[1, 1], xlabel = L"X_1", ylabel = L"X_2", limits = (-8, 0, -5, 5))

    lines!(ax, [-8.0, 0.0], [-2.0, 6.0], color = :black, linestyle = :dash, label = L"X_1 - X_2 = -6")

    scatter!(ax, X_conditional[1, :], X_conditional[2, :], color = :red, markersize = 7, label = L"X \sim \mathcal{N}(\mu, \Sigma) \mid (X_1 - X_2 < -6)")
    scatter!(ax, X_unconditional[1, :],  X_unconditional[2, :],  color = :blue, markersize = 7, label = L"X \sim \mathcal{N}(\mu, \Sigma)")

    axislegend(ax, position = :lt)
    fig
end

References

1. Z. I. Botev. The normal law under linear restrictions: Simulation and estimation via minimax tilting. Journal of the Royal Statistical Society. Series B, Statistical methodology 79, 125–148 (2017).