Constrained Inference

As discussed in the main manuscript^[1], we show that the adaptive algorithm (AGESS) is ergodic for fairly general target distributions on $(\mathcal{X}, \mathcal{B}(\mathcal{X}))$, where $\mathcal{X}$ is an open subset of $\mathbb{R}^P$. This allows us to consider using AGESS in constrained inference settings that can posed as inference on an open set. Consider the following constrained inference linear regression problem^[2]:

\[p(\boldsymbol{\beta} \mid \mathbf{y}) \propto \exp(-0.5 \lVert \mathbf{y} - \mathbf{X}\boldsymbol{\beta}\rVert^2_2) \mathbb{1}_{\{\boldsymbol{\beta}: \lVert \boldsymbol{\beta} \rVert^2_2 < 1\}}(\boldsymbol{\beta}).\]

Thus, we can see that the support is constrained to the interior of the unit ball – an open set. We will start by generating data from an unconstrained linear regression framework, with $\boldsymbol{\beta}$ outside the set of interest.

import Random
using AdaptEllipticalSliceSampler
using Distributions
using Plots
using LinearAlgebra

Random.seed!(123)

D = 2
N = 100

### Generate Design matrix
X = randn(N, D)
### True beta outside of set of interest
β = [-1, 1]
### Generate Response Variable
y = X * β .+ randn(N)

100-element Vector{Float64}:
  0.1025289544904367
  3.365201258498111
 -0.1606258787847763
  0.7015626928584558
  0.6874634999812321
 -0.6996806134788154
 -0.36518887820634693
  2.19779727005064
 -0.24081661546943667
 -1.162239327931908
  ⋮
 -2.2392257213354245
 -3.0318512872643724
  0.23620724223442757
 -0.7072244137885617
 -0.7603916657170077
 -2.067320130521191
 -4.348747698330303
  4.208679046545524
 -1.5206513106811579

Now that we have generated our data, let's construct a function that evaluates the posterior log pdf (up to a constant). Note: This function can be implemented more efficiently, however this implementation is sufficient for this demonstration.

function log_posterior(β::AbstractVector{Y}, X::AbstractMatrix{Y},
                       y::AbstractVector{Y}) where {Y <: AbstractFloat}
    l_lpdf::Float64 = -Inf
    ### Incorporate constraint
    if norm(β) < 1
        l_lpdf = -0.5 * norm(y - X * β)^2
    end

    return l_lpdf
end

log_posterior (generic function with 1 method)

With this function, we can use the AGESS function to draw samples from the posterior.

### Specify the dimension of the target distribution
P = D
### Specify the number of MCMC iterations
n_MCMC = 5000

### Run AGESS
results = AGESS(β -> log_posterior(β, X, y), n_MCMC, P)

AdaptEllipticalSliceSampler.MCMC_output([0.0 0.0; -0.0818748226318659 0.04027972475090179; … ; -0.7864407265100615 0.5527953645451524; -0.6239835598194547 0.7431858312026026], [-126.43695591242208, -117.19507486598643, -112.9913601603847, -106.60893184847751, -102.08389035230314, -102.89148407621472, -102.17565290079962, -98.21945640930792, -98.40691470391232, -96.61634922574157  …  -63.494957586885455, -63.9986789447422, -63.00301217984904, -63.722923357570366, -63.41227690424116, -63.08250044857275, -62.47855477261998, -62.79060995941325, -63.054732962559385, -62.705369003349354], AGESS_MCMC_params(Main.var"#2#3"(), [0.0, 0.0], [1.0 0.0; 0.0 1.0], true, 6.0, 0.5, 0.05, 0.5, 0.05, 5000, 2), [0.0042059678213805605 0.0031061039209897383; 0.0031061039209897383 0.004527979560298203], [-0.7053042201950414, 0.6483631759691028])

Let's visualize the samples from the target distribution.

function circleShape()
    θ = LinRange(0, 2*π, 100)
    x = cos.(θ)
    y = sin.(θ)
    return x, y
end

plot(circleShape(), ylims=(-1.5, 1.5), xlims=(-1.5, 1.5), legend=false)

scatter!([β[1]], [β[2]], color="red")
scatter!(results.samps[1000:end,1], results.samps[1000:end,2], alpha = 0.4)

In the figure above, we can see that the red point is the true unconstrained $\boldsymbol{\beta}$, from which the data was generated from. In green, we can see samples generated from our Markov chain. We can see that AGESS provides an alternative method to conduct constrained Bayesian inference – targeting the exact target distribution, rather than smooth relaxations of the sharply constrained priors ^{[2] [3] [4]}.