HistSmoother

Documentation for the histogram smoother of Wand and Yu (2022).

This model is available through the BayesDensityHistSmoother package.

The HistSmoother approach is based on modelling the data-generating density as a log-spline with a fixed and relatively large basis dimension $K$, with a smoothness-inducing prior distribution on the spline coefficients. More specifically, the model specification is given by the hierarchical scheme

\[\begin{align*} x_i\,|\, \boldsymbol{\beta} &\sim \frac{\exp\{\sum_{k=1}^K \beta_k z_k(x_i)\}}{\int \exp\{\sum_{k=1}^K \beta_k z_k(x)\}\,\text{d}x}, &i = 1,\ldots, n,\\ \beta_k &\sim \mathrm{Normal}(0, \sigma_\beta^2), &k = 1, 2,\\ \beta_k \, |\, \sigma^2 &\sim \mathrm{Normal}(0, \sigma^2), &k \geq 3\\ \sigma^2 &\sim \text{HalfCauchy}(0, s_\sigma), \end{align*}\]

where $z_k$ are O'Sullivan spline basis functions (Wand and Ormerod, 2008) and $\sigma_\beta, s_\beta > 0$ are fixed hyperparameters.

The estimates are then obtained by smoothing a likelihood based on the Poisson approximation to a histogram.

This module implements the Gibbs sampling approach of Wand and Yu (2022) for Markov chain Monte Carlo-based inference. We also provide an implementation of the semiparametric mean-field variational inference algorithm proposed in the same paper.

Module API

BayesDensityHistSmoother.HistSmoother — Type

HistSmoother{T<:Real} <: AbstractBayesDensityModel{T}

Struct representing a spline histogram smoother model.

Constructors

HistSmoother(x::AbstractVector{<:Real}; kwargs...)
HistSmoother{T}(x::AbstractVector{<:Real}; kwargs...)

Arguments

x: The data vector.

Keyword arguments

K: B-spline basis dimension of a regular augmented spline basis. Defaults to 52.
bounds: A tuple giving the support of the B-spline mixture model.
n_bins: Number of bins used to construct the histogram likelihood. Defaults to 400.
prior_scale_fixed: Scale hyperparameter for 0th and 1st order (fixed effect) spline terms. Defaults to 1000.0.
prior_scale_random: Scale hyperparameter for the (higher order) mixed effects spline terms. Defaults to 1000.0.

Returns

shs: A histogram spline smoother object.

Examples

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

julia> shs = HistSmoother(x)
52-dimensional HistSmoother{Float64}:
Using 5000 binned observations with 400 bins.
 support: (-0.105, 1.105)
Hyperparameters:
 prior_scale_fixed = 1000.0
 prior_scale_random = 1000.0

julia> shs = HistSmoother(x; K = 80, prior_scale_fixed = 1e5);

Extended help

Binning

The binning step used by the spline histogram smoother is an essential part of the model fitting procedure, and can as such not be disabled. Using a greater number of bins means that less precision is lost due to the binning step, but makes the model fitting procedure slower due to a larger compuatational burden. Note that the number of bins only affects the model fitting process, and does otherwise not change the returned

Hyperparameter selection

The hyperparameter s_β contols the smoothness of the resulting density estimates. Setting this to a smaller value leads to smoother estimates.

source

Evaluating the pdf and cdf

Distributions.pdf — Method

pdf(
    shs::HistSmoother,
    params::NamedTuple,
    t::Union{Real, AbstractVector{<:Real}}
) -> Union{Real, Vector{<:Real}}

pdf(
    shs::HistSmoother,
    params::AbstractVector{NamedTuple},
    t::Union{Real, AbstractVector{<:Real}}
) -> Matrix{<:Real}

Evaluate $f(t\, |\, \boldsymbol{\eta})$ for a given HistSmoother when the model parameters of the NamedTuple params are given by $\boldsymbol{\eta}$.

The named tuple should contain a field named :β. If the parameters argument does not contain a field named :norm, then the normalization constant will be computed using Simpson's method. Alternatively, if parameters contains the field :norm, then this value is used instead.

source

Distributions.cdf — Method

cdf(
    shs::HistSmoother,
    params::NamedTuple,
    t::Union{Real, AbstractVector{<:Real}}
) -> Union{Real, Vector{<:Real}}

cdf(
    shs::HistSmoother,
    params::AbstractVector{NamedTuple},
    t::Union{Real, AbstractVector{<:Real}}
) -> Matrix{<:Real}

Evaluate $F(t\, |\, \boldsymbol{\eta})$ for a given HistSmoother when the model parameters of the NamedTuple params are given by $\boldsymbol{\eta}$.

Internally, this function computes the cdf on a predefined regular grid, and uses linear interpolation to approximate the cdf.

source

Utility functions

BayesDensityCore.hyperparams — Method

hyperparams(
    shs::HistSmoother{T}
) where {T} -> @NamedTuple{prior_scale_fixed::T, prior_scale_random::T}

Returns the hyperparameters of the spline histogram smoother shs as a NamedTuple.

source

Markov chain Monte Carlo

StatsBase.sample — Method

sample(
    [rng::Random.AbstractRNG],
    hs::HistSmoother{T},
    n_samples::Int;
    n_burnin::Int              = min(100, div(n_samples, 5)),
    initial_params::NamedTuple = get_default_initparams_mcmc(hs)
) where {T} -> PosteriorSamples{T}

Generate n_samples posterior samples from a HistSmoother using an augmented Gibbs sampler.

Arguments

rng: Optional random seed used for random variate generation.
hs: The HistSmoother object for which posterior samples are generated.
n_samples: The total number of samples (including burn-in).

Keyword arguments

n_burnin: Number of burn-in samples.
initial_params: Initial values used in the MCMC algorithm. Should be supplied as a NamedTuple with fields :β and :σ2, where :β is a K-dimensional vector and σ2 is a positive scalar.

Returns

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

Examples

julia> using Random

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

julia> hs = HistSmoother(x);

julia> ps = sample(Xoshiro(1), hs, 1100);

source

Variational inference

BayesDensityHistSmoother.HistSmootherVIPosterior — Type

HistSmootherVIPosterior{T<:Real} <: AbstractVIPosterior{T}

Struct representing the variational posterior distribution of a HistSmoother.

Fields

q_β: Distribution representing the optimal variational density q*(β).
q_σ: Distribution representing the optimal variational density q*(σ²).
shs: The HistSmoother to which the variational posterior was fit.

source

BayesDensityCore.varinf — Method

varinf(
    hs::HistSmoother{T};
    initial_params::NamedTuple = _get_default_initparams_varinf(hs),
    max_iter::Int              = 500,
    rtol::Real                 = 1e-5
) where {T} -> HistSmootherVIPosterior{T}

Find a variational approximation to the posterior distribution of a HistSmoother using semiparametric mean-field variational inference.

Arguments

hs: The HistSmoother whose posterior we want to approximate.

Keyword arguments

initial_params: Initial values of the VI parameters μ_opt, Σ_opt and b_σ_opt.
max_iter: Maximal number of VI iterations. Defaults to 500.
rtol: Relative tolerance used to determine convergence. Defaults to 1e-5.

Returns

vip: A HistSmootherVIPosterior object representing the variational posterior.
info: A VariationalOptimizationResult describing the result of the optimization.

Note

To run the optimization loop for a fixed number of iterations irrespective of the convergence criterion, one can set rtol = 0.0, and max_iter equal to the desired total iteration count. Note that setting rtol to a strictly negative value will issue a warning.

Examples

julia> using Random

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

julia> hs = HistSmoother(x);

julia> vip, info = varinf(hs; rtol=1e-6);

Extended help

Convergence

The criterion used to determine convergence is that the relative change in the expectation of $\mathbb{E}(\sigma^{-2})$ falls below the given rtol.

source