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