Introduction¶

The Bayesian Quantile Regression (BQReg) library is a lightweight C++ implementation of Bayesian quantile regression using the asymmetric Laplace density representation of the problem, with bindings for Python and R.

Features:

A header-only C++11 library with OpenMP-accelerated MCMC sampling for parallel computation.
Built on the Eigen (version >= 3.4.0) templated linear algebra library for fast and efficient matrix-based computation.
Straightforward linking with parallelized BLAS libraries, such as OpenBLAS.
Seamless integration with Python and R through Pybind11 and Rcpp modules.
Available as a single precision (float) or double precision (double) library.
Released under a permissive, non-GPL license.

Author: Keith O’Hara

License: Apache Version 2.0

Installation¶

The C++ implementation of BQReg is a header-only library. First, clone the library and related submodules:

# clone optim into the current directory
git clone https://github.com/kthohr/BayesianQuantileRegression ./BayesianQuantileRegression

# change directory
cd ./BayesianQuantileRegression

# clone necessary submodules
git submodule update --init

Then simply add the BQReg header files (found under cpp/include) to your project using:

#include "bqreg.hpp"

Beyond the files contained in the Git submodules, the only dependencies are a copy of Eigen (version >= 3.4.0) and a C++14-compatible compiler.

For detailed instructions on how to install the R and Python bindings, see the installation page.

Contents¶

Installation¶

The C++ implementation of BQReg is a header-only library, and we provide bindings for Python and R.

C++¶

The C++ library requires only the Eigen(3) C++ linear algebra library, as well as the submodules included in the source code. (Note that Eigen version 3.4.0 requires a C++14-compatible compiler.)

First clone the library and required submodules:

# clone optim into the current directory
git clone https://github.com/kthohr/BayesianQuantileRegression ./BayesianQuantileRegression

# change directory
cd ./BayesianQuantileRegression

# clone necessary submodules
git submodule update --init

Then simply add the header files (found under cpp/include) to your project using:

#include "bqreg.hpp"

The following options should be declared before including the BQReg header files.

OpenMP functionality is enabled by default if the _OPENMP macro is detected (e.g., by invoking -fopenmp with GCC or Clang).
- To explicitly enable OpenMP features, use:
```
#define BQREG_USE_OPENMP
```
- To explicitly disable OpenMP functionality, use:
```
#define BQREG_DONT_USE_OPENMP
```

Python¶

If you do not already have Eigen3 and Pybind11, clone and set environment variables:

# clone Eigen3
git clone https://gitlab.com/libeigen/eigen.git

# clone Pybind11
git clone https://github.com/pybind/pybind11.git

# set environment variables
export EIGEN_INCLUDE_PATH = "<path to this directory>/eigen"
export PYBIND11_INCLUDE_PATH = "<path to this directory>/pybind11/include"

Then install from the Python subdirectory:

# change directory into the Python subdirectory
cd BayesianQuantileRegression/Python

# install the package
python3 -m pip install . --user

R¶

First install the RcppEigen package from R:

install.packages("RcppEigen")

Then install the package from the R subdirectory:

# change directory into the R subdirectory
cd BayesianQuantileRegression/R

# install the package
R CMD INSTALL .

Model Description¶

Classical Quantile Regression¶

Define the quantile function of a random variable, \(Y\), conditional on a random vector, \(X\), as:

\[Q_{Y | X} (\tau) := \inf \left\{ y : F_{Y|X}(y) \geq \tau \right\}\]

where \(\tau \in [0,1]\) and \(F_{Y|X}\) denotes the conditional distribution function of \(Y\) given \(X\).

Quantile regression models the conditional quantile function of a continuous random variable as an affine function of a feature set, \(X\), and a parameter vector, \(\beta\). That is: \(Q_{Y | X} (\tau) = X \beta(\tau)\), where we emphasize that the parameter vector, \(\beta\), varies by the quantile of interest.

The standard frequentist estimator for \(\beta(\tau)\) is given by the solution to the following optimization problem:

\[\hat{\beta}(\tau) = \arg \min_{\beta} \left\{ \sum_{i=1}^N \rho_{\tau} (Y_i - X_i \beta) \right\}\]

where \(\rho_{\tau}\) denotes the check function:

\[\rho_{\tau} (y) = y (\tau - \mathbf{1}[y < 0])\]

In practice, this optimization problem is solved by linear programming techniques.

Bayesian Quantile Regression¶

The Bayesian version of the classical quantile regression problem assumes that the quantile “error” follows an Asymmetric Laplace Distribution (ALD), the density function of which is given by:

\[f(y | \mu, \sigma, \tau) = \frac{\tau (1 - \tau)}{\sigma} \exp \left( - \rho_{\tau} \left( \frac{y - \mu}{\sigma} \right) \right)\]

where \(\mu\) is the location parameter, \(\sigma\) is the scale parameter, and \(\tau\) is the “asymmetry” parameter. Note that when \(\tau = 0.5\), the ALD density collapses to that of the Laplace distribution (with a rescaled scaling parameter).

The statistical model is defined using a mean-variance mixture representation of the ALD:

\[\begin{aligned} Y_i &= X_i \beta (\tau) + \theta(\tau) \sigma z_i + \omega(\tau) \sigma \sqrt{z_i} u_i \\ u_i &\stackrel{\text{iid}}{\sim} N(0,1) \\ z_i &\stackrel{\text{iid}}{\sim} E(1) \end{aligned}\]

where the quantile-dependent parameters \(\theta(\tau)\) and \(\omega(\tau)\) are defined as:

\[\begin{aligned} \theta(\tau) &= \frac{1 - 2\tau}{\tau (1 - \tau)} \\ \omega(\tau) &= \sqrt{ \frac{2}{\tau (1 - \tau)} } \end{aligned}\]

The likelihood function is then given by:

\[L(\beta, \sigma | Y, X, Z, \tau) \propto \prod_{i=1}^N \frac{1}{\omega \sigma \sqrt{z_i}} \exp \left( - \frac{(Y_i - X_i \beta - \theta \sigma z_i )^2}{2 \omega^2 \sigma^2 z_i} \right)\]

The Gibbs sampling procedure of Kozumi and Kobayashi (2011) proceeds in three blocks.

Draw \(\beta\):

\[\begin{aligned} \beta | Y, X, Z, \sigma &\sim N(\widetilde{\beta}, \widetilde{V}) \\ \widetilde{\beta} &= \widetilde{V} \left( V_0^{-1} \beta_0 + \frac{1}{\omega^2 \sigma^2} \sum_{i=1}^N \frac{1}{z_i} X_i^\top (Y_i - \theta \sigma z_i) \right) \\ \widetilde{V}^{-1} &= V_0^{-1} + \frac{1}{\omega^2 \sigma^2} \sum_{i=1}^N \frac{1}{z_i} X_i^\top X_i \end{aligned}\]

Draw \(\{ z_i \}_{i=1}^N\):

\[\begin{aligned} z_i | Y, X, \beta, \sigma &\sim g(x; a, b) \propto \frac{1}{\sqrt{x}} \exp \left( - \frac{1}{2} (a^2 x + b_i^2 x^{-1}) \right) \\ a &= \sqrt{ \frac{2}{\sigma^2} + \frac{\theta^2}{\sigma^2 \omega^2} } \\ b_i &= \frac{|Y_i - X_i \beta|}{\sigma \omega} \end{aligned}\]

The sampling density for \(z_i\) is proportional to a generalized inverse-Gaussian distribution.

Draw \(\sigma\):

\[\begin{aligned} \sigma | Y, X, Z, \beta &\sim IG(\tilde{\gamma}_0, \tilde{\gamma}_1) \\ \tilde{\gamma}_0 &= \gamma_0 + \frac{3}{2} N \\ \tilde{\gamma}_1 &= \gamma_1 + \sum_{i=1}^N \left[ \sigma z_i + \frac{1}{2 \omega^2 \sigma z_i} ( Y_i - X_i \beta - \theta \sigma z_i ) \right] \end{aligned}\]

C++¶

Class Declaration¶

class bqreg_t¶

Bayesian Quantile Regression class.

Public Functions

~bqreg_t() = default¶: Default destructor.

bqreg_t(const bqreg_t &obj_inp)¶

Constructor using a bqreg_t object to copy from.

Parameters: obj_inp – an object of type bqreg_t whose elements will be used to initalize a new object of the same type.

bqreg_t(bqreg_t &&obj_inp)¶

Constructor using a moveable bqreg_t object.

Parameters: obj_inp – a moveable object of type bqreg_t whose elements will be used to initalize a new object of the same type.

inline explicit bqreg_t(const ColVec_t &Y_inp, const Mat_t &X_inp)¶

Constructor using data matrices.

Parameters

Y_inp – an n x 1 vector defining the target variable
X_inp – an n x K matrix of features

inline explicit bqreg_t(ColVec_t &&Y_inp, Mat_t &&X_inp)¶

Constructor using moveable data matrices.

Parameters

Y_inp – an n x 1 vector defining the target variable
X_inp – an n x K matrix of features

inline bqreg_t &operator=(const bqreg_t &obj_inp)¶

Assignment operator.

Parameters: obj_inp – an object of type bqreg_t whose elements will be used to initalize a new object of the same type.
Returns: an object of type bqreg_t

inline bqreg_t &operator=(bqreg_t &&obj_inp)¶

Assignment operator.

Parameters: obj_inp – a moveable object of type bqreg_t whose elements will be used to initalize a new object of the same type.
Returns: an object of type bqreg_t

inline int get_omp_n_threads()¶

OpenMP threads check.

Returns: the number of OpenMP threads to use in the Gibbs sampler stage.

inline void set_omp_n_threads(const int omp_n_threads_inp)¶

OpenMP threads set.

Parameters: omp_n_threads_inp – the number of OpenMP threads to use in the Gibbs sampler stage.

inline void set_seed_value(const size_t seed_val_inp)¶

RNG engine seeding.

Parameters: seed_val_inp – the seed value for the RNG engine.

inline void load_data(const ColVec_t &Y_inp, const Mat_t &X_inp)¶

Load data.

Parameters

Y_inp – an n x 1 vector defining the target variable
X_inp – an n x K matrix of features

inline void set_quantile_target(const fp_t tau_inp)¶

Set the target quantile value.

Parameters: tau_inp – a real value between zero and one

inline void set_prior_params(const ColVec_t &prior_beta_mean_inp, const Mat_t &prior_beta_var_inp, const fp_t prior_sigma_shape_inp, const fp_t prior_sigma_scale_inp)¶

Set the Prior.

Set the parameters of the distributions that define the prior

Parameters

prior_beta_mean_inp – mean of the prior distribution for \( \beta \)
prior_beta_var_inp – variance of the prior distribution for \( \beta \)
prior_sigma_shape_inp – shape parameter of the prior distribution for \( \sigma \)
prior_sigma_scale_inp – scale parameter of the prior distribution for \( \sigma \)

inline ColVec_t get_initial_beta_draw()¶

Initial value check for the Gibbs sampler.

Returns: a vector defining the initial values for \( \beta \) in the Gibbs sampler.

inline void set_initial_beta_draw(const ColVec_t &beta_initial_draw_inp)¶

Set the initial values for the Gibbs sampler.

Parameters: beta_initial_draw_inp – a vector defining the initial values for \( \beta \) in the Gibbs sampler.

inline void gibbs(const size_t n_burnin_draws, const size_t n_keep_draws, const size_t thinning_factor, Mat_t &beta_draws, Mat_t &z_draws, ColVec_t &sigma_draws)¶

Run the Gibbs sampler.

Parameters

n_burnin_draws – the number of burnin draws
n_keep_draws – the number of draws to keep, post burnin
thinning_factor – the number of draws to skip between keep draws
beta_draws – a writable matrix to store the draws of \( \beta \)
z_draws – a writable matrix to store the draws of \( z \)
sigma_draws – a writable vector to store the draws of \( \sigma \)

Public Members

ColVec_t Y¶: An n x 1 vector defining the target variable

Mat_t X¶: An n x K matrix of features

fp_t tau¶: The target quantile value

ColVec_t prior_beta_mean¶: Mean of the prior distribution for \( \beta \)

Mat_t prior_beta_var¶: Variance of the prior distribution for \( \beta \)

fp_t prior_sigma_shape¶: Shape parameter of the prior distribution for \( \sigma \)

fp_t prior_sigma_scale¶: Scale parameter of the prior distribution for \( \sigma \)

Example¶

Code to run a median regression example is given below.

#include "bqreg.hpp"

inline
Eigen::MatrixXd
eigen_randn(const size_t nr, const size_t nc)
{
    static std::mt19937 gen{ std::random_device{}() };
    static std::normal_distribution<> dist;

    return Eigen::MatrixXd{ nr, nc }.unaryExpr([&](double x) { (void)(x); return dist(gen); });
}

int main()
{
    // generate data
    const int n = 500;
    const int K = 3;

    Eigen::MatrixXd X = eigen_randn(n, K);
    X.col(0).setConstant(1);

    Eigen::VectorXd beta0(3);
    beta0 << 5.0, 1.3, 1.8;

    Eigen::VectorXd Y = (X * beta0).array() + eigen_randn(n, 1).array();

    // initialize object
    bqreg::bqreg_t obj = bqreg::bqreg_t(Y, X);

    // set prior pars
    Eigen::VectorXd beta_bar = Eigen::VectorXd::Zero(K);

    Eigen::MatrixXd V0 = Eigen::MatrixXd::Zero(K, K);
    V0.diagonal().array() = 1.0 / 0.001;

    int prior_shape = 3;
    int prior_scale = 3;

    obj.set_prior_params(beta_bar, V0, prior_shape, prior_scale);

    // (optional) set the number of OpenMP threads
    obj.set_omp_n_threads(4);
    std::cout << "Number of OpenMP threads: " << obj.get_omp_n_threads() << std::endl;

    // (optional) set the seed value of the Gibbs sampler RNG
    obj.set_seed_value(1111);

    // set the target quantile
    obj.set_quantile_target(0.5);

    // run Gibbs sampler
    size_t n_burnin_draws = 10000;
    size_t n_keep_draws = 10000;
    size_t thinning_factor = 0;

    Eigen::MatrixXd beta_draws;
    Eigen::MatrixXd z_draws;
    Eigen::VectorXd sigma_draws;

    obj.gibbs(n_burnin_draws, n_keep_draws, thinning_factor, beta_draws, z_draws, sigma_draws);

    //

    std::cout << "\n - beta posterior mean:\n" << beta_draws.rowwise().mean() << std::endl;

    //

    return 0;
}

Python¶

Class Declaration¶

class BayesianQuantileRegression¶

Bayesian Quantile Regression class

Public Functions

__init__(self, Union[pd.Series, pd.DataFrame, np.ndarray] target, Union[pd.Series, pd.DataFrame, np.ndarray] features)¶

Initialize the BayesianQuantileRegression class

    Parameters:
        target: An n x 1 vector defining the target variable (Y)
        features: An n x K matrix of features (X)

set_seed_value(self, int seed_value)¶

Set the RNG seed value for Gibbs sampling

    Parameters:
        seed_value: An integer value

get_omp_n_threads(self)¶

Get the current value determining the number of OpenMP threads to use with the Gibbs sampler

    Returns:
        An integer value

set_omp_n_threads(self, int omp_n_threads)¶

Set the number of OpenMP threads to use for the Gibbs sampler

    Parameters:
        omp_n_threads: An integer value

set_prior_params(self, np.ndarray beta_mean, np.ndarray beta_var, float sigma_shape, float sigma_scale)¶

Set the parameters of the prior distributions

    Parameters:
        beta_mean: Mean vector of the normal prior for beta
        beta_var: Variance matrix of the normal prior for beta
        sigma_shape: Shape parameter of the prior distribution for sigma^2
        sigma_scale: Scale parameter of the prior distribution for sigma^2

get_initial_beta_draw(self)¶

Return the initial values of beta for the Gibbs sampler

    Parameters:
        beta_initial_draw: Initial draw vector

set_initial_beta_draw(self, np.ndarray beta_initial_draw)¶

Set the initial values of beta for the Gibbs sampler

    Parameters:
        beta_initial_draw: Initial draw vector

fit(self, float tau=0.5, int n_burnin_draws=1000, int n_keep_draws=1000, int thinning_factor=0)¶

Fit method for the BayesianQuantileRegression class

    Parameters:
        tau: the target quantile value
        n_burnin_draws: the number of burn-in draws
        n_keep_draws: the number of post burn-in draws to return
        thinning_factor: the number of draws to skip between keep draws
    
    Returns:
        A tuple of matrices containing posterior draws, ordered as follows: (beta, z, sigma)
    
    Notes:
        The total number of draws will be: n_burnin_draws + (thinning_factor + 1) * n_keep_draws

Example¶

Code to run a median regression example is given below.

# import libraries
import numpy as np
import matplotlib.pyplot as plt
from pybqreg import BayesianQuantileRegression

# generate data
n = 500
K = 3

X = np.random.normal(0.0, 1.0, [n, K])
X[:,0] = np.ones(n)

beta0 = np.random.uniform(1.0, 2.0, K)
beta0[0] = 5.0

Y = np.matmul(X, beta0) + np.random.normal(0.0, 1.0, n)

# initialize object
obj = BayesianQuantileRegression(Y,X)

# set prior pars
beta_bar = np.zeros(K)

V0 = np.zeros([K, K])
np.fill_diagonal(V0, 1.0/0.001)

prior_shape = 3
prior_scale = 3

obj.set_prior_params(beta_bar, V0, prior_shape, prior_scale)

# (optional) set the initial draw for beta
beta_hat = np.linalg.solve( np.matmul(X.transpose(),X), np.matmul(X.transpose(),Y) )
obj.set_initial_beta_draw(beta_hat)

# (optional) set the number of OpenMP threads
obj.set_omp_n_threads(4)
obj.get_omp_n_threads()

# (optional) set the RNG seed value of the Gibbs sampler
obj.set_seed_value(1111)

# set the target quantile (tau) and run the Gibbs sampler
tau = 0.5

n_burnin_draws = 10000
n_keep_draws = 10000
thinning_factor = 0

beta_draws, z_draws, sigma_draws = obj.fit(tau, n_burnin_draws, n_keep_draws, thinning_factor)

beta_mean = np.mean(beta_draws, axis = 1)
beta_std = np.std(beta_draws, axis = 1)

# np.mean(z_draws, axis = 1) / np.mean(sigma_draws)

# plot beta draws

fig, axs = plt.subplots(1, 3, figsize=(15, 5), tight_layout=True)

for k in range(3):
    axs[k].hist(beta_draws[k,:], bins = 200)

    axs[k].set_xlim(beta_mean[k] - 3*beta_std[k], beta_mean[k] + 3*beta_std[k])

    axs[k].axvline(beta_mean[k], c = 'r')

plt.show()

R¶

Class Declaration¶

TBW

Example¶

Code to run a median regression example is given below.

library(BQReg.Rcpp)

set.seed(1900)

#

n <- 1000
K <- 3

X <- matrix(rnorm(n*K),n,K)
X[,1] <- 1

beta0 <- matrix(runif(K,1,2), ncol = 1)
beta0[1] <- 5.0

Y <- X %*% beta0 + rnorm(n)

#

bqreg_obj <- new(bqreg)

bqreg_obj$load_data(Y, X)

#

beta_bar <- matrix(0, K, 1)

Vbar <- diag(K) * 1/0.001

n0 <- 3
s0 <- 3

bqreg_obj$set_prior_params(beta_bar, Vbar, n0, s0)

#

beta_hat <- solve( t(X) %*% X, t(X) %*% Y )
bqreg_obj$set_initial_beta_draw(beta_hat)

#

tau <- 0.5

bqreg_obj$set_quantile_target(tau)

n_burnin_draws <- 10000
n_keep_draws <- 10000

gibbs_res <- bqreg_obj$gibbs(n_burnin_draws, n_keep_draws, 0)

rowMeans(gibbs_res$beta_draws)