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()