General API

This page explains how to fit the Bayesian density models implemented in BayesDensity.jl. Most of the methods implemented in this package support two modes of posterior inference: simulation consistent inference through Markov chain Monte Carlo (MCMC) and approximate through variational inference (VI). We also document most of the convenience methods available for computing select posterior quantities of interest, such as the posterior mean or quantiles of $f(t)$ for some $t \in \mathbb{R}$.

The plotting API of this package is documented on a separate page

Defining models

The first step to estimating a density with this package is to create a model object for which posterior inference is desired. All density models in this package are subtypes of AbstractBayesDensityModel:

BayesDensityCore.AbstractBayesDensityModel — Type

AbstractBayesDensityModel{T<:Real}

Abstract super type for all Bayesian density models implemented in this package.

In order to create a model object, we call the corresponding contructor with the data and other positional- and keyword arguments. For example, we can create a BSplineMixture object with default hyperparameters as follows:

bsm = BSplineMixture(randn(1000))

For more detailed information on the arguments supported by each specific Bayesian density model we refer the reader to the methods documentation.

Evaluating the density and the cumulative distribution function

The density estimators implemented in this package all specify a model $f(t\,|\, \boldsymbol{\eta})$ for the density of the data, which depends on a parameter $\boldsymbol{\eta}$. In order to calculate $f(\cdot)$ for a given $\boldsymbol{\eta}$, each Bayesian density model implements the pdf method.

Distributions.pdf — Method

pdf(
    bdm::AbstractBayesDensityModel,
    parameters::NamedTuple,
    t::Union{Real, AbstractVector{<:Real}}
) -> Union{Real, Vector{<:Real}}

pdf(
    bdm::AbstractBayesDensityModel,
    parameters::AbstractVector{<:NamedTuple},
    t::Union{Real, AbstractVector{<:Real}}
) -> Matrix{<:Real}

Evaluate $f(t \,|\,\boldsymbol{\eta})$ of the Bayesian density model bdm for every $\boldsymbol{\eta}$ in parameters and every element in the collection t.

If a single NamedTuple is passed to the parameters argument, this function outputs either a scalar or a vector depending on the input type of the third argument t.

If a vector of NamedTuples is passed to the second positional argument, then this function returns a Matrix of size (length(t), length(parameters)).

For models that only implement the signature pdf(::AbstractBayesDensityModel, ::Any, ::Real), a generic fallback method is provided for vectors of parameters and vector evaluation grids. However, it is recommended that most models provide specialized methods for vectors of parameters and vectors of evaluation points, as it is often possible to implement batch evaluation more efficiently (e.g. by leveraging BLAS calls instead of loops) when the parameters and the evaluation grid are provided in batches.

The cumulative distribution function of a model can be computed in a similar way by using the cdf method:

Distributions.cdf — Method

cdf(
    bdm::AbstractBayesDensityModel,
    parameters::NamedTuple,
    t::Union{Real, AbstractVector{<:Real}}
) -> Union{Real, Vector{<:Real}}

cdf(
    bdm::AbstractBayesDensityModel,
    parameters::AbstractVector{<:NamedTuple},
    t::Union{Real, AbstractVector{<:Real}}
) -> Matrix{<:Real}

Evaluate the cumulative distribution function $F(t\,|\, \boldsymbol{\eta}) = \int_{-\infty}^t f(s\,|\,\boldsymbol{\eta})\,\text{d}s$ of the Bayesian density model bdm for every $\boldsymbol{\eta}$ in parameters and every element in the collection t.

If a single NamedTuple is passed to the parameters argument, this function outputs either a scalar or a vector depending on the input type of the third argument t.

If a vector of NamedTuples is passed to the second positional argument, then this function returns a Matrix of size (length(t), length(parameters)).

Generic fallback methods for computing the cdf for vectors of parameters and vector evaluation grids are also provided for models that implement the signature cdf(::AbstractBayesDensityModel, ::Any, ::Real).

Other methods

All of the density models implemented in this package depend on the choice of various hyperparameters, which can be retrieved by utilizing the following method:

BayesDensityCore.hyperparams — Method

hyperparams(bdm::AbstractBayesDensityModel) -> @NamedTuple

Return the hyperparameters of the model bdm as a NamedTuple.

For the exact format of the returned hyperparameters for a specific Bayesian density model type, we refer to the docs of the individual density estimators.

To compute the support of a given model, the support method is provided.

Distributions.support — Method

support(bdm::AbstractBayesDensityModel{T}) where {T} -> NTuple{2, T}

Return the support of the model bdm as an 2-dimensional tuple.

The element type of the model object can be determined via the eltype method:

Base.eltype — Method

eltype(::AbstractBayesDensityModel{T}) where {T}

Return the element type of a Bayesian density model.

eltype(::PosteriorSamples{T}) where {T}

Get the element type of a PosteriorSamples object.

eltype(::AbstractVIPosterior{T}) where {T}

Get the element type of a variational posterior object.

The following function is used to select a grid used for plotting of fitted model objects. A default fallback is provided, but it may be necessary to overload this method when implementing new models in order to make plotting functions work without having to supply an explicit grid of values for the first coordinate axis.

BayesDensityCore.default_grid_points — Function

default_grid_points(bdm::AbstractBayesDensityModel{T}) where {T} -> AbstractVector{T}

Get the default grid used for plotting of density estimates.

Defaults to returning constructing a grid based on the extrema of bdm.data.x. If a given struct does not store a copy of the original data used to construct the model object as bdm.data.x, this method should be implemented.

Markov chain Monte Carlo

The main workhorse of MCMC-based inference is the sample method, which takes a Bayesian density model object as input and generates posterior samples through a specialized MCMC routine.

StatsBase.sample — Method

sample(
    [rng::Random.AbstractRNG],
    bdm::AbstractBayesDensityModel{T},
    n_samples::Int
    args...;
    [n_burnin::Int],
    kwargs...
) where {T} -> PosteriorSamples{T}

Generate approximate posterior samples from the density model bdm using Markov chain Monte Carlo methods.

This functions returns a PosteriorSamples object which can be used to compute posterior quantities of interest such as the posterior mean of $f(t)$ or posterior quantiles. See the specific method docstring for method-specific additional positional and keyword arguments.

Arguments

rng: Seed used for random variate generation.
bdm: The Bayesian density model object to generate posterior samples from.
n_samples: Number of Monte Carlo samples (including burn-in)
args...: Other model-specific positional arguments. See the specific model docstring for this method for further details.

Keyword arguments

n_burnin: Number of burn-in samples.
kwargs...: Other model-specific kwargs. See the specific model docstring for this method for further details.

Returns

ps: A PosteriorSamples object holding the posterior samples and the original model object.

All of the implemented MCMC methods return an object of type PosteriorSamples:

BayesDensityCore.PosteriorSamples — Type

PosteriorSamples{T<:Real}

Struct holding posterior samples of the parameters of a Bayesian density model.

Fields

samples: Vector holding posterior samples of model parameters.
model: The model object to which samples were fit.
n_samples: Total number of Monte Carlo samples.
non_burnin_ind: Indices of non-burnin samples.

BayesDensityCore.samples — Method

samples(ps::PosteriorSamples{T, V}) where {T, V<:AbstractVector} -> V

Get the posterior samples of a PosteriorSamples object as a vector.

The following methods can be used to extract useful information about the model object, such as the underlying model object and the number of samples.

BayesDensityCore.model — Method

model(ps::PosteriorSamples) -> AbstractBayesDensityModel

Return the model object of ps.

BayesDensityCore.n_samples — Method

n_samples(ps::PosteriorSamples) -> Int

Get the total number of samples of a PosteriorSamples object, including burn-in samples.

BayesDensityCore.n_burnin — Method

n_burnin(ps::PosteriorSamples) -> Int

Get the number of burn-in samples of a PosteriorSamples object.

By default, PosteriorSamples objects also store the burn-in samples from the MCMC routine. These can be discarded via the following method:

BayesDensityCore.drop_burnin — Method

drop_burnin(ps::PosteriorSamples{T}) where {T} -> PosteriorSamples{T}

Create a new PosteriorSamples object where the burn-in samples have been discarded.

Multiple PosteriorSamples objects can also be concatenated to create a single PosteriorSamples object. This is particularly useful when a preliminary MCMC run is deemed to be too short, and one wants to pool the original samples with the samples from a new MCMC run.

Base.vcat — Method

vcat(ps::PosteriorSamples...) -> PosteriorSamples

Concatenate the samples of multiple PosteriorSamples objects to form a single object.

Computing posterior summary statistics

When using Bayesian density estimators, we are often interested in computing various summary statistics of the posterior draws from an MCMC procedure. For instance, we may be interested in providing an estimate of the density $f$ (e.g. the posterior mean) and to quantify the uncertainty in this estimate (e.g. via credible bands).

To this end, BayesDensityCore provides methods for PosteriorSamples objects that let us easily compute relevant summary statistics for the density $f$, as shown in the short example below:

bsm = BSplineMixture(randn(1000))
posterior = sample(bsm, 2000; n_burnin=400)

# Compute the posterior mean of f(0.5)
mean(posterior, pdf, 0.5)

# Compute the posterior 0.05 and 0.95-quantiles of f(0.5)
# Note that supplying pdf as the second argument is optional here
quantile(posterior, 0.5, [0.05, 0.95]) == quantile(posterior, pdf, 0.5, [0.05, 0.95])

In some cases it may also be of interest to carry out posterior inference for the cumulative distribution function $F(t) = \int_{-\infty}^t f(s)\, \text{d}s$. Computing posterior summary statistics for the cdf instead of the pdf is easily achieved by replacing the pdf in the second argument with cdf instead:

# Compute the posterior mean of F(0.5)
mean(posterior, cdf, 0.5)

# Compute the posterior 0.05 and 0.95-quantiles of F(0.5)
# Note that supplying cdf as the second argument is necessary here
quantile(posterior, cdf, 0.5, [0.05, 0.95])

The posterior summary statistics available through BayesDensityCore are the following:

Statistics.mean — Method

mean(
    ps::PosteriorSamples,
    [func = pdf],
    t::Union{Real, AbstractVector{<:Real}}
) -> Union{Real, Vector{<:Real}}

Compute the approximate posterior mean of a functional of $f$ for every element in the collection t using Monte Carlo samples.

The target functional can be either be the pdf $f$ or the cdf $F$, and is controlled by adjusting the func argument. By default, the posterior mean of $f$ is computed.

If the input t is a scalar, a scalar is returned. If t is a vector, this function returns a vector the same length as t.

Examples

julia> x = (1.0 .- (1.0 .- LinRange(0, 1, 5001)) .^(1/3)).^(1/3);

julia> ps = sample(Random.Xoshiro(1), BSplineMixture(x), 5000);

julia> mean(ps, 0.1)
0.0969450407517681

julia> mean(ps, [0.1, 0.8])
2-element Vector{Float64}:
 0.0969450407517681
 1.3662358915400654

Statistics.quantile — Method

quantile(
    ps::PosteriorSamples,
    [func = pdf],
    t::Union{Real, AbstractVector{<:Real}},
    q::RealAbstractVector{<:Real,
) -> Union{Real, Vector{<:Real}}

quantile(
    ps::PosteriorSamples,
    [func = pdf]
    t::Union{Real, AbstractVector{<:Real}},
    q::AbstractVector{<:Real},
) -> Matrix{<:Real}

Compute the approximate posterior quantile(s) of a functional of $f$ for every element in the collection t using Monte Carlo samples.

The target functional can be either be the pdf $f$ or the cdf $F$, and is controlled by adjusting the func argument. By default, the posterior quantiles of $f$ are computed.

In the case where both t and q are scalars, the output is a real number. When t is a vector and q a scalar, this function returns a vector of the same length as t. If q is supplied as a Vector, then this method returns a Matrix of dimension (length(t), length(q)), where each column corresponds to a given quantile. This is also the case when t is supplied as a scalar.

Examples

julia> x = (1.0 .- (1.0 .- LinRange(0, 1, 5001)) .^(1/3)).^(1/3);

julia> ps = sample(Random.Xoshiro(1), BSplineMixture(x), 5000);

julia> quantile(ps, 0.1, 0.5); # Get the posterior median of f(0.5)

julia> quantile(ps, cdf, 0.1, 0.5); # Get the posterior median of F(0.5)

julia> quantile(ps, [0.1, 0.8], [0.05, 0.95]); # Get the posterior 0.05, 0.95-quantiles of f(0.1) and f(0.8)

Statistics.median — Method

median(
    ps::PosteriorSamples,
    [func = pdf]
    t::Union{Real, AbstractVector{<:Real}},
) -> Union{Real, Vector{<:Real}}

Compute the approximate posterior median of a functional of $f$ for every element in the collection t using Monte Carlo samples.

The target functional can be either be the pdf $f$ or the cdf $F$, and is controlled by adjusting the func argument. By default, the posterior median of $f$ is computed.

If the input t is a scalar, a scalar is returned. If t is a vector, this function returns a vector the same length as t.

Statistics.var — Method

var(
    ps::PosteriorSamples,
    [func = pdf],
    t::Union{Real, AbstractVector{<:Real}}
) -> Union{Real, Vector{<:Real}}

Compute the approximate posterior variance of a functional of $f$ for every element in the collection t using Monte Carlo samples.

The target functional can be either be the pdf $f$ or the cdf $F$, and is controlled by adjusting the func argument. By default, the posterior variance of $f$ is computed.

If the input t is a scalar, a scalar is returned. If t is a vector, this function returns a vector the same length as t.

Examples

julia> x = (1.0 .- (1.0 .- LinRange(0, 1, 5001)) .^(1/3)).^(1/3);

julia> ps = sample(Random.Xoshiro(1), BSplineMixture(x), 5000);

julia> var(ps, 0.1) # get the posterior variance of f(0.1)
0.00027756364767372627

julia> var(ps, [0.1, 0.8]) # get the posterior variance of f(0.1) and f(0.8)
2-element Vector{Float64}:
 0.00027756364767372627
 0.005977674286240125

Statistics.std — Method

std(
    ps::PosteriorSamples,
    [func = pdf],
    t::Union{Real, AbstractVector{<:Real}}
) -> Union{Real, Vector{<:Real}}

Compute the approximate posterior standard deviation of a functional of $f$ for every element in the collection t using Monte Carlo samples.

The target functional can be either be the pdf $f$ or the cdf $F$, and is controlled by adjusting the func argument. By default, the posterior standard deviation of $f$ is computed.

If the input t is a scalar, a scalar is returned. If t is a vector, this function returns a vector the same length as t. This method is equivalent to sqrt.(var(rng, vip, t, n_samples)).

It is also possible to evaluate the pdf or cdf for all non-burn in samples at a grid of input points:

Distributions.pdf — Method

pdf(
    ps::PosteriorSamples,
    t::Union{<:Real, AbstractVector{<:Real}}
) -> Matrix{<:Real}

Evaluate the pdf $f(t \,|\,\boldsymbol{\eta})$ of the Bayesian density model to which ps was fit for every sample $\boldsymbol{\eta}^{(s)}$ from the posterior and every element in the collection t.

This function always returns a Matrix of size (length(t), n_samples(ps) - n_burnin(ps)).

Distributions.cdf — Method

cdf(
    ps::PosteriorSamples,
    t::Union{<:Real, AbstractVector{<:Real}}
) -> Matrix{<:Real}

Evaluate the cdf $F(t \,|\,\boldsymbol{\eta})$ of the Bayesian density model to which ps was fit for every sample $\boldsymbol{\eta}^{(s)}$ from the posterior and every element in the collection t.

This function always returns a Matrix of size (length(t), n_samples(ps) - n_burnin(ps)).

Variational inference

The varinf method can be used to compute a variational approximation to the posterior distribution:

BayesDensityCore.varinf — Method

varinf(
    bdm::AbstractBayesDensityModel{T},
    args...;
    kwargs...
) where {T}

Compute a variational approximation to the posterior distribution.

The positional arguments and keyword arguments supported by this function, as well as the type of the returned variational posterior object differs between different subtypes of AbstractBayesDensityModel.

Any call to varinf will return a subtype of the abstract type AbstractVIPosterior:

BayesDensityCore.AbstractVIPosterior — Type

AbstractVIPosterior{T<:Real}

Abstract super type representing the variational posterior distribution of AbstractBayesDensityModel

For most models, varinf also returns an object which stores the result of the optimization procedure, see VariationalOptimizationResult.

The following convenience methods are also part of the public API:

BayesDensityCore.model — Method

model(vip::AbstractVIPosterior{T}) where {T} -> AbstractBayesDensityModel{T}

Get the model object to which the variational posterior vip was fitted.

Base.eltype — Method

eltype(::AbstractBayesDensityModel{T}) where {T}

Return the element type of a Bayesian density model.

eltype(::PosteriorSamples{T}) where {T}

Get the element type of a PosteriorSamples object.

eltype(::AbstractVIPosterior{T}) where {T}

Get the element type of a variational posterior object.

Generating samples from the variational posterior

The sample method makes it possible to generate independent samples from the variational posterior. This is particularly useful in cases where inference for multiple posterior quantities (e.g. medians, variances) is desired.

StatsBase.sample — Method

sample(
    [rng::Random.AbstractRNG],
    vip::AbstractVIPosterior{T},
    n_samples::Int
) where {T} -> PosteriorSamples{T}

Generate n_samples independent samples from the variationonal posterior distribution vip.

Examples

julia> x = (1.0 .- (1.0 .- LinRange(0, 1, 5001)) .^(1/3)).^(1/3);

julia> vip = varinf(BSplineMixture(x));

julia> vps = sample(Random.Xoshiro(1812), vip, 5000);

As shown in the above docstring, using the sample method on a AbstractVIPosterior object returns an object of type PosteriorSamples. As such, all of the convenience methods showcased in the previous subsection will also work for the object returned by sample.

Computing posterior summary statistics

BayesDensityCore also provides convenience methods for AbstractVIPosterior objects that let us easily compute relevant summary statistics for the density $f$ and the cdf $F$ directly from the variational posterior object:

bsm = BSplineMixture(randn(1000))
viposterior, info = varinf(bsm)

# Compute the (variational) posterior mean of f(0.5)
mean(viposterior, pdf, 0.5)

# Compute the (variational) posterior median of F(0.5)
median(viposterior, cdf, 0.5)

# Compute the (variational) posterior 0.05 and 0.95-quantiles of f(0.5)
# Note that supplying pdf as the second argument is optional here
quantile(viposterior, 0.5, [0.05, 0.95]) ≈ quantile(viposterior, pdf, 0.5, [0.05, 0.95])

The full list of available summary statistics is the same as that for PosteriorSamples objects:

Statistics.mean — Method

mean(
    [rng::Random.AbstractRNG],
    vip::AbstractVIPosterior,
    [func = pdf],
    t::Union{Real, AbstractVector{<:Real}},
    [n_samples::Int=1000]
) -> Union{Real, Vector{<:Real}}

Compute the approximate posterior mean of $f(t)$ for every element in the collection t using Monte Carlo samples.

If the input t is a scalar, a scalar is returned. If t is a vector, this function returns a vector the same length as t.

Examples

julia> x = (1.0 .- (1.0 .- LinRange(0, 1, 5001)) .^(1/3)).^(1/3);

julia> vip = varinf(BSplineMixture(x));

julia> mean(Random.Xoshiro(1), vip, 0.1)
0.08615412808594237

julia> mean(Random.Xoshiro(1), vip, [0.1, 0.8])
2-element Vector{Float64}:
 0.08615412808594237
 1.3674886390998342

Statistics.quantile — Method

quantile(
    [rng::Random.AbstractRNG],
    vip::AbstractVIPosterior,
    [func = pdf],
    t::Union{Real, AbstractVector{<:Real}},
    q::Union{Real, AbstractVector{<:Real}},
    [n_samples::Int=1000]
) -> Union{Real, Vector{<:Real}}

quantile(
    [rng::Random.AbstractRNG],
    vip::AbstractVIPosterior,
    [func = pdf],
    t::Union{Real, AbstractVector{<:Real}},
    q::AbstractVector{<:Real},
    [n_samples::Int=1000]
) -> Matrix{<:Real}

Compute the posterior q-quantile(s) of the pdf $f(t)$ or the cdf $F(t)$ for each element in the collection t.

By default this function falls back to quantile(sample(rng, vip, n_samples), func, t, q)

In the case where both t and q are scalars, the output is a real number. When t is a vector and q a scalar, this function returns a vector of the same length as t. If q is supplied as a Vector, then this function returns a Matrix of dimension (length(t), length(q)), where each column corresponds to a given quantile. This is also the case when t is supplied as a scalar.

Examples

julia> x = (1.0 .- (1.0 .- LinRange(0, 1, 5001)) .^(1/3)).^(1/3);

julia> vip = varinf(BSplineMixture(x));

julia> quantile(Random.Xoshiro(1), vip, 0.9, 0.5)
0.537450082172813

julia> quantile(Random.Xoshiro(1), vip, [0.2, 0.8], [0.05, 0.95])
2×2 Matrix{Float64}:
 0.34042  0.362478
 1.3039   1.43599

Statistics.median — Method

median(
    [rng::Random.AbstractRNG],
    vip::AbstractVIPosterior,
    [func = pdf],
    t::Union{Real, AbstractVector{<:Real}},
    [n_samples::Int=1000]
) -> Union{Real, Vector{<:Real}}

Compute the posterior median of the pdf $f(t)$ or the cdf $F(t)$ for each element in the collection t.

Equivalent to quantile(rng, vip, t, 0.5, n_samples).

In the case where both t and q are scalars, the output is a real number. When t is a vector, this function returns a vector of the same length as t.

Statistics.var — Method

var(
    [rng::Random.AbstractRNG],
    vip::AbstractVIPosterior,
    [func = pdf],
    t::Union{Real, AbstractVector{<:Real}},
    [n_samples::Int=1000]
) -> Union{Real, Vector{<:Real}}

Compute the posterior variance of the pdf $f(t)$ or the cdf $F(t)$ for each element in the collection t.

If the input t is a scalar, a scalar is returned. If t is a vector, this function returns a vector the same length as t.

Examples

julia> x = (1.0 .- (1.0 .- LinRange(0, 1, 5001)) .^(1/3)).^(1/3);

julia> vip = varinf(BSplineMixture(x));

julia> var(Random.Xoshiro(1), vip, 0.1)
3.2895012264929507e-6

julia> var(Random.Xoshiro(1), vip, [0.1, 0.8])
2-element Vector{Float64}:
 3.2895012264929507e-6
 0.0014793006688108977

Statistics.std — Method

std(
    [rng::Random.AbstractRNG],
    vip::AbstractVIPosterior,
    [func = pdf],
    t::Union{Real, AbstractVector{<:Real}},
    [n_samples::Int=1000]
) -> Union{Real, Vector{<:Real}}

Compute the posterior standard deviation of the pdf $f(t)$ or the cdf $F(t)$ for each element in the collection t.

If the input t is a scalar, a scalar is returned. If t is a vector, this function returns a vector the same length as t. This method is equivalent to sqrt.(var(rng, vip, t, n_samples)).

Note

Note that each call to mean, quantile, median, var or std in most cases will first simulate a random sample from the posterior distribution, and then uses this sample to compute a Monte Carlo approximation of the quantity of interest using these samples. If posterior inference for multiple quantities is desired, then it is recommended to first use sample, and call these functions on this object as only a single batch of posterior samples is generated in this case.

Storing info from the variational optimization

In order to provide a simple way of performing convergence diagnostics for variational optimization problems, BayesDensityCore exports the VariationalOptimizationResult type.

BayesDensityCore.VariationalOptimizationResult — Type

VariationalOptimizationResult{T<:Real}

Struct holding the result of a variational inference procedure.

Fields

ELBO: The values of the evidence lower bound per iteration.
converged: Boolean flag indicating whether the optimization was succesful or not.
n_iter: Number of iterations run before termination.
tolerance: Tolerance parameter used to determine convergence.
variational_posterior: The fitted variational posterior distribution.

BayesDensityCore.elbo — Function

elbo(varoptinf::VariationalOptimizationResult{T}) where {T} -> Vector{T}

Get the value of the evidence lower bound for each iteration of the optimization procedure.

BayesDensityCore.n_iter — Function

n_iter(varoptinf::VariationalOptimizationResult{T}) where {T} -> Int

Get the number of iterations used to fit the variational posterior distribution.

BayesDensityCore.converged — Function

converged(varoptinf::VariationalOptimizationResult{T}) where {T} -> Bool

Return the convergence status of the variational optimization as a bool.

BayesDensityCore.tolerance — Function

tolerance(varoptinf::VariationalOptimizationResult{T}) where {T} -> T

Get the value of the tolerance level used to determine convergence.

BayesDensityCore.posterior — Function

posterior(varoptinf::VariationalOptimizationResult{T}) where {T} -> AbstractVIPosterior{T}

Get the fitted variational posterior distribution.