Bayesian Linear Regression in Julia
Author willtebbutt
13 Stars
Updated Last
1 Year Ago
Started In
March 2019

Bayesian Linear Regression in Julia

Stable Dev Build Status Codecov

This is a simple package that does one thing, Bayesian Linear Regression, in around 100 lines of code.

Intended Use and Functionality

The interface sits at roughly the same level as that of Distributions.jl. This means that while you won't find a scikit-learn-style fit function, you will find all of the primitives that you need to construct such a function to suit your particular problem. In particular, one can:

  • Construct a BayesianLinearRegressor (BLR) object by providing a mean-vector and precision matrix for the weights of said regressor. This object represents a distribution over (linear) functions.
  • "Index into" said distribution over functions to construct an IndexedBLR object, which represents a finite-dimensional marginal of a BayesianLinearRegressor.
  • Compute the log marginal likelihood of a vector of observations.
  • Sample from the finite-dimensional marginals of a BLR.
  • Perform posterior inference to produce a new BLR with an update mean and precision.
  • All operations are fully compatible with Zygote.jl (hopefully), so you can use gradient-based optimisation to tune the hyperparameters of your regressor etc.

For examples of how to use this package in conjunction with Flux and Zygote, see the examples directory.


A BayesianLinearRegressor in D dimensions works with data where:

  • inputs X should be a D x N matrix of Reals where each column is from one data point.
  • outputs y should be an N-vector of Reals, where each element is from one data point.

Example Usage

# Install the packages if you don't already have them installed
] add BayesianLinearRegressors LinearAlgebra Random Optim Plots Zygote
using BayesianLinearRegressors, LinearAlgebra, Random, Optim, Plots, Zygote

# Fix seed for re-producibility.
rng = MersenneTwister(123456)

# We don't export anything, so you need to explicitly import the stuff that you need.
using BayesianLinearRegressors: BayesianLinearRegressor, logpdf, rand, posterior, marginals,

# Construct a BayesianLinearRegressor prior over linear functions of `X`.
mw, Λw = zeros(2), Diagonal(ones(2))
f = BayesianLinearRegressor(mw, Λw)

# Index into the regressor and assume heterscedastic observation noise `Σ_noise`.
N = 10
X = collect(hcat(collect(range(-5.0, 5.0, length=N)), ones(N))')
Σ_noise = Diagonal(exp.(randn(N)))
fX = f(X, Σ_noise)

# Generate some toy data by sampling from the prior.
y = rand(rng, fX)

# Compute the adjoint of `rand` w.r.t. everything given random sensitivities of y′.
_, back_rand = Zygote.pullback(
    (X, Σ_noise, mw, Λw)->rand(rng, BayesianLinearRegressor(mw, Λw)(X, Σ_noise), 5),
    X, Σ_noise, mw, Λw,
back_rand(randn(N, 5))

# Compute the `logpdf`. Read as `the log probability of observing `y` at `X` under `f`, and
# Gaussian observation noise with zero-mean and covariance `Σ_noise`.
logpdf(fX, y)

# Compute the gradient of the `logpdf` w.r.t. everything.
    (X, Σ_noise, y, mw, Λw)->logpdf(BayesianLinearRegressor(mw, Λw)(X, Σ_noise), y),
    X, Σ_noise, y, mw, Λw,

# Perform posterior inference. Note that `f′` has the same type as `f`.
f′ = posterior(fX, y)

# Compute `logpdf` of the observations under the posterior predictive (because why not?)
logpdf(f′(X, Σ_noise), y)

# Sample from the posterior predictive distribution.
N_plt = 1000
X_plt = hcat(collect(range(-6.0, 6.0, length=N_plt)), ones(N_plt))'
f′X_plt = rand(rng, f′(X_plt, eps()), 100) # Samples with machine-epsilon noise for stability

# Compute some posterior marginal statisics.
normals = marginals(f′(X_plt, eps()))
m′X_plt = mean.(normals)
σ′X_plt = std.(normals)

# Plot the posterior marginals.
plotly(); # My prefered backend. Use a different one if you prefer / this doesn't work.
posterior_plot = plot();
plot!(posterior_plot, X_plt[1, :], f′X_plt; # Posterior samples.
plot!(posterior_plot, X_plt[1, :], [m′X_plt m′X_plt]; # Posterior credible intervals.
    fillrange=[m′X_plt .- 3 .* σ′X_plt, m′X_plt .+ 3 * σ′X_plt],
plot!(posterior_plot, X_plt[1, :], m′X_plt; # Posterior mean.
scatter!(posterior_plot, X[1, :], y; # Observations.

Up For Grabs

  • Scikit-learn style interface: it wouldn't be too hard to implement a scikit-learn - style interface to handle basic regression tasks, so please feel free to make a PR that implements this.
  • Monte Carlo VI (MCVI): i.e. variational inference using the reparametrisation trick. This could be very useful when working with large data sets and applying big non-linear transformations, such as neural networks, to the inputs as it would enable mini-batching. I would envise at least supporting both a dense approximate posterior covariance and diagonal (i.e. mean-field), where the latter is for small-moderate dimensionalities and the latter for very high-dimensional problems.

Bugs, Issues, and PRs

Please do report and bugs you find by raising an issue. Please also feel free to raise PRs, especially if for one of the above Up For Grabs items. Raise an issue to discuss the extension in detail before opening a PR if you prefer though.

Related Work

BayesianLinearRegression.jl is closely related, but appears to be a WIP and hasn't been touched in around a year or so (as of 27-03-2019).

Used By Packages

No packages found.